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ABSTRACT 

We present Star Formation Histories (SFHs) for a sample of 104 massive (stellar mass 
M > 10 1() Mg) quiescent galaxies (MQGs) at z — 1.0-1.5 from the analysis of spectro- 
photometric data from the SHARDS and HST/WFC3 G102 and G141 surveys of the GOODS- 
N field, jointly with broad-band observations from ultraviolet (UV) to far-infrared (Far-IR). 
The sample is constructed on the basis of rest-frame UVJ colours and specific star forma¬ 
tion rates (sSFR=SFR/Mass). The Spectral Energy Distributions (SEDs) of each galaxy are 
compared to models assuming a delayed exponentially declining SFH. A Monte Carlo al¬ 
gorithm characterizes the degeneracies, which we are able to break taking advantage of the 
SHARDS data resolution, by measuring indices such as MgUV and D4000. The population 
of MQGs shows a duality in their properties. The sample is dominated (85%) by galaxies with 
young mass-weighted ages, 1m <2 Gyr, short star formation timescales, (r) ~ 60-200 Myr, 
and masses log(M/M 0 ) ~ 10.5. There is an older population (15%) with 1 m-2 -4 Gyr, longer 
star formation timescales, (r) ~ 400 Myr, and larger masses, log(M/M Q ) ~ 10.7. The SFHs 
of our MQGs are consistent with the slope and the location of the Main Sequence (MS) of 
star-forming galaxies at z > 1.0, when our galaxies were 0.5-1.0 Gyr old. According to these 
SFHs, all the MQGs experienced a Fuminous Infrared Galaxy (FIRG) phase that lasts for 
~ 500 Myr, and half of them an Ultra Fuminous Infrared Galaxy (UFIRG) phase for ~ 100 
Myr. We find that the MQG population is almost assembled at z ~ 1, and continues evolving 
passively with few additions to the population. 

Key words: galaxies: formation; galaxies: evolution; galaxies: high-redshift; galaxies: stellar 
content 


1 INTRODUCTION 

In the current paradigm of cosmic evolution, galaxies grow their 
mass by accreting gas from the cosmic web (e.g. Tacconi et al. 
2010) and transforming it into stars. A tight relation between mass 
* E-mail: helenads@ucm.es and SFR exists for normal star-forming galaxies, known as the MS 
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(e.g. Noeske et al. 2007; Elbaz et al. 2007; Rodighiero et al. 2010). 
Galaxies grow in mass within or above this MS until, eventually, the 
exhaustion of gas or some feedback mechanism halt the star for¬ 
mation. At this time, galaxies reach a quiescence state. The mass 
evolution of the quiescent population is then dominated by (dry) 
mergers of already assembled (smaller or similar in mass) galaxies. 
Lambda Cold Dark Matter (ACDM) theory requires more massive 
halos to generally assemble later than less massive ones. However, 
over the past few decades there have been a series of observational 
results indicating that many processes, such as star formation, oc¬ 
cur earlier in the most massive and luminous galaxies than in less 
massive galaxies, in a scenario called downsizing (e.g., Cowie et al. 
1996; Madau et al. 1996; Cimatti et al. 2006; Bundy et al. 2007; 
Perez-Gonzalez et al. 2008; Ilbert et al. 2010). 

Focusing on the evolution of the most massive galax¬ 
ies along the lifetime of the Universe, we now know that 
many of such systems were already assembled at high redshifts 
(Perez-Gonzalez et al. 2008; Marchesini et al. 2009; Ilbert et al. 
2013; Muzzin et al. 2013; Tomczaketal. 2014; Grazian et al. 
2015). Moreover, a significant fraction of high-z massive galax¬ 
ies were not actively forming stars and were evolving pas¬ 
sively (Cimatti et al. 2004; Daddi et al. 2005; Papovich et al. 2006; 
Fontana et al. 2009; Santini et al. 2009; Dominguez Sanchez et al. 
2011; Ilbert et al. 2013), i.e., they were MQGs. The number den¬ 
sities of MQGs at intermediate redshifts are in disagreement 
with semi-analytical models (e.g. Pozzetti et al. 2010; Ilbert et al. 
2013; Muzzin et al. 2013). Besides, the observed MQGs at high- 
z are found to be much more compact than their local analogues, 
implying a strong mass-size evolution with cosmic time (e.g. 
Trujillo et al. 2006; Buitrago et al. 2008; van Dokkum et al. 2008; 
Barro et al. 2014). 

Understanding the formation mechanisms of this population 
of MQGs is fundamental to improve our picture of galaxy evo¬ 
lution. In particular, having an accurate description of the SFHs 
of these galaxies is crucial to have a good estimation of the time 
needed by physical processes to ignite the star formation and then 
quench it in massive galaxies. However, up to date, there are few 
results on the individual properties and SFHs of MQGs at high-z, 
mostly because estimating SFHs is very hard with the data we typ¬ 
ically have at hand. 

Quiescent galaxies at high-z are difficult to observe in the op¬ 
tical bands as their emission in the rest-frame UV is weak (due 
to the absence of star formation processes), and dominated by 
absorption features (e.g., Daddi et al. 2005; Cimatti et al. 2008; 
van Dokkum et al. 2011). Indeed, these galaxies show very weak 
or no emission lines (Kriek et al. 2009; Belli et al. 2014) and ab¬ 
sorption features related to an old and passively evolved stellar 
population. Spectroscopic observations are time consuming and are 
therefore limited to a small number of galaxies or to the analysis of 
stacked spectra (Kriek et al. 2009; Onodera et al. 2012; Toft et al. 
2012; Bedregal et al. 2013; Whitaker et al. 2013; Onodera et al. 
2015; Belli et al. 2015). 

Thanks to the arrival of deep and wide multi-wavelength pho¬ 
tometric surveys such as COSMOS (Scoville et al. 2007), ULTRA¬ 
VISTA (McCracken et al. 2012) or CANDELS (Grogin et al. 2011; 
Koekemoer et al. 2011), the mass functions and number densities 
of quiescent galaxies have been studied up to z ~ 3-4. But the spec¬ 
tral resolution of photometric data is not enough to study in detail 
the stellar population properties. The main reason is the presence of 
strong degeneracies in the analysis of their SEDs using stellar pop¬ 
ulation synthesis models. These degeneracies are mainly related to 
the similar effect of different levels of dust attenuation, age, and 


metallicity in very low spectral resolution data (R ~ 7, typical of 
broad-band studies). This indeed complicates the accurate determi¬ 
nation of high-z galaxy properties based on SED-fitting; being the 
stellar mass the only reliable parameter obtained in these kind of 
studies (Eisner et al. 2008; Santini et al. 2015). 

In order to break the typical degeneracies inherent to any study 
of stellar populations in distant galaxies, we need data with higher 
spectral resolution than broad-band photometry (Pacifici et al. 
2013). There are spectral features which help to break these degen¬ 
eracies. For example, the Mgyy index probes several absorption 
lines (e.g., Mgl T2852, Mgll TT2796, 2804, Fell TT2600, 2606) 
and has been shown to be extremely reliable to identify galaxies 
dominated by evolved stars. Moreover, these absorption lines can 
be used to easily distinguish the SED of a MQG from the feature¬ 
less spectrum of a dusty starburst (Daddi et al. 2005). The break 
in the stellar continuum at 4000 A, D4000, is also a good age in¬ 
dicator (Bruzual 1983; Balogh et al. 1999; Kauffmann et al. 2003). 
It arises because of the accumulation of a large number of spec¬ 
tral lines in a narrow wavelength region. Given the dependence of 
D4000 on stellar atmospheric parameters (Gorgas et al. 1999), it is 
very prominent for galaxies older than ~ 1 Gyr. The strength of 
the 4000 A break for a single stellar population increases with its 
age and depends weakly on the metallicity at low ages « 1 Gyr, 
Hernan-Caballero et al. 2013). The Mguv and D4000 indices have 
been successfully used in the past to obtain redshifts and ages 
of stellar populations in massive galaxies at high-z (Saracco et al. 
2005; Kriek et al. 2011; Ferreras et al. 2012). 

In this work, we take advantage of the spectro-photometric 
resolution of the Survey for High-z Absorption Red and Dead 
Sources survey (SHARDS, Perez-Gonzalez et al. 2013) to obtain 
robust estimations of the SFH of MQGs at high-z. SHARDS is 
an ultra-deep optical survey of the GOODS-N field covering the 
wavelength range between 500 and 950 nm with 25 contiguous 
medium-band filters, providing a spectral resolution R ~ 50. We 
combine the SHARDS observations with spectroscopic data from 
the HST WFC3 G102 and G141 grisms covering 900-1600 nm, as 
well as with multi-wavelength ancillary data from UV to Far-IR 
extracted from the Rainbow database (Perez-Gonzalez et al. 2008; 
Barro etal. 2011). We, therefore, use data with a spectral resolu¬ 
tion R ~ 50 or better from 500 nm to 1600 nm (jointly with broad¬ 
band photometry), a significant improvement from previous works 
on the subject. We perform a SED-fitting to the whole wavelength 
range (up to the IRAC bands) using delayed exponentially declin¬ 
ing stellar population models. Thanks to the unique photometric 
dataset and spectral coverage of this work, we are able to account 
for the degeneracies and study in detail and individually the SFHs 
(masses, SFRs, ages, star formation timescales, dust attenuations 
and metallicities) of a sample of MQGs at z =1.0 - 1.5, discussing 
the implications for the early mass assembly of galaxies. 

The structure of this paper is as follows: in Sect. 2, we present 
the data and the sample selection. In Sect. 3, we explain the SED- 
fitting method and the procedure to characterize the degeneracies 
inherent to the stellar population analysis. In Sect. 4, we analyse the 
derived properties and the time-evolution of our sample of MQGs 
on the basis of their SFHs. Finally, in Sect. 5, we summarize our 
conclusions. 

Throughout this work, we standardize to a (h, Cl M , fL\)=(0.7, 
0.3, 0.7) Wilkinson Microwave Anisotropy Probe (WMAP) 
concordance cosmology (Spergel et al. 2003), AB magnitudes, 
(Oke & Gunn 1983), a Kroupa 2001 IMF (integrated from 0.1 - 
100 Mo) and Bruzual & Chariot (2003) (BC03, hereafter) stellar 
population synthesis models. 
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2 DATA AND SAMPLE SELECTION 

In this work, we analyse a sample of MQGs (M > 10 10 M 0 ) at 1.0 < 
z < 1.5 selected with two different criteria: their rest-frame UVJ 
colours and their sSFR. In the following subsections, we describe 
the datasets gathered for this work, as well as the details about the 
selection criteria. 


2.1 Data 

Our sample was selected in the 130 arcmin 2 covered by SHARDS 
(Perez-Gonzalez et al. 2013) in the GOODS-N region. SHARDS is 
an ESO/GTC Large Program carried out with the OSIRIS instru¬ 
ment on the 10.4 m Gran Telescopio Canarias (GTC). It consists of 
an ultra-deep optical spectro-photometric survey of the GOODS-N 
field at wavelengths between 500 and 950 nm using 25 contiguous 
medium-band filters which provide a spectral resolution R ~ 50. 
The data reach an AB magnitude of 26.5 (at least at a 3cr level) 
with sub-arcsec seeing in each one of the 25 bands. More details 
about the reduction and calibration of the SHARDS data can be 
found in Perez-Gonzalez et al. (2013). 

For this paper, we complement the SHARDS data with the 
G102 and G141 grism observations of the GOODS-N field car¬ 
ried out with the HST/WFC3 instrument. The G102 spectroscopic 
program (PI: Barro) covers the spectral region between 800 and 
1150 nm with R ~ 210 and a 5cr magnitude limit of 21.5 mag in 
the F140W band. The G141 data (PI: Weiner) covers from 1100 to 
1700 nm with R ~ 130 and a 5cr limit of 21.5 mag in the F160W 
band. The WFC3 grism spectroscopic data were reduced using the 
aXe software version 2.3. Based on the 2D spectra provided by 
aXe, we extracted ID spectra for each galaxy in our sample using 
the effective radius as the extraction width. We visually inspected 
each spectrum, adjusting the extraction parameters (aperture width 
and spectral range), to avoid contamination from nearby sources. 
We refer the reader to Esquej et al. (in preparation) for a detailed 
description of the reduction and extraction of WFC3 grism data. 

lointly with the SHARDS and WFC3 grism data, we also 
use in our analysis the ancillary multi-wavelength catalogue 
and advanced products in the GOODS-N field presented in 
Perez-Gonzalez et al. (2008) and compiled in the Rainbow Cos¬ 
mological Surveys Database 1 (see also Perez-Gonzalez et al. 2008; 
Barro et al. 2011). This dataset includes observations from X-rays 
to the Far-IR and radio bands, as well as spectroscopic data in 
the GOODS-N field from the literature. In Perez-Gonzalez et al. 
(2008), merged photometric catalogues were presented, including 
broad-band data for a stellar mass selected sample based on ultra¬ 
deep IRAC observations. For this work, we have merged this cata¬ 
logue with the SHARDS and WFC3 grism datasets. Using the full 
SED, we carried out a stellar population modeling which provided 
accurate photometric redshifts, stellar masses, SFRs, and rest- 
frame synthetic colours for ~ 26,000 stellar mass selected galaxies 
in the region of the GOODS-North field covered by SHARDS. 

Concerning photometric redshifts, the ultra-deep medium 
band SHARDS data allowed us to obtain high quality photometric 
redshifts for all sources, which will be presented in Barro et al. (in 
preparation). The median |Az|/(l+z) is 0.0067 for the 2650 sources 
with I < 25 and spectroscopic redshifts (see also Ferreras et al. 
2014). 


1 http://rainbowx.fis.ucm.es 


As described in Perez-Gonzalez et al. (2008), the parent sam¬ 
ple used in this paper, selected with ultra-deep IRAC imaging as a 
proxy for stellar mass, is complete for galaxies with M > 10 9 5 M 0 
up to z = 1.5 and a maximally old stellar population. The IRAC 
selection is biased against younger galaxies with large attenuations 
and masses M < 10 10 M 0 , so we impose this mass limit in the defi¬ 
nition of our final sample of MQGs. 

Synthetic rest-frame colours were estimated for all galaxies 
in the parent sample by convolving the best-fitting stellar popula¬ 
tion models with transmission curves for standard filters. In this 
paper, we will use the lohnson U and V filters, as well as the 
2MASS 7-band filter to construct a UVJ diagram and select quies¬ 
cent galaxies. We note that the actual transmission curves for which 
we estimated UVJ rest-frame absolute magnitudes were taken from 
EAZY (Brammer et al. 2008) and FAST (Kriek et al. 2009) 2 , in 
order to match the colour distribution in the UVJ diagrams used 
in Whitaker et al. (2011), from which we extracted the quiescence 
definition. 

Finally, in order to select quiescent galaxies for this paper, 
SFRs were estimated for all galaxies in the parent sample based 
on either mid- and far-IR data from Spitzer and Herschel, or from 
UV luminosities. In the case of IR emitters, GOODS-N has been 
observed with the deepest MIPS data in the sky, with a 5 cr limit of 
30 pjy (Perez-Gonzalez et al. 2005), and also very deep PACS and 
SPIRE images (Elbaz et al. 2011), with the following 5 cr limits: 
1.7 and 3.6 mly for PACS 100 and 160 pm, and 9, 12, and 11 mjy 
for SPIRE 250, 350, and 500 pm (these limits include the effect of 
confusion in the SPIRE bands). 

We looked for counterparts for the galaxies in our mass se¬ 
lected sample in the Spitzer MIPS bands using a search radius of 
2" for 24 pm, as described in Perez-Gonzalez et al. (2005, 2008). 
For Herschel PACS and SPIRE bands, we built merged photometric 
catalogues using position priors from the MIPS data, as described 
in Perez-Gonzalez et al. (2010) and Rawle et al. (2015). The Her¬ 
schel catalogues are then linked to the MIPS catalog, so we used 
the same search radius. We only considered detections above the 
5tr flux limit for each band. 

Infrared Luminosities, h IR , were estimated from the Spitzer 
and Herschel data by fitting dust emission models from 
Chary & Elbaz (2001) to all available photometric data points 
with rest-frame wavelength longer than 6 pm. We also checked 
that similar results (with typical differences smaller than 0.2 dex) 
were obtained using the templates from Rieke et al. (2009) and 
Dale & Helou (2002). When we only had a single photometric 
data point (in almost all cases, 24 /im) we scaled the models from 
Chary & Elbaz (2001) to the monochromatic luminosity probed by 
that observation. We compared our IR-based SFRs with those ob¬ 
tained by applying the method described in Wuyts et al. (2011); 
Rujopakarn et al. (2012) and Elbaz et al. (2011), obtaining a good 
(< 0.2 dex) agreement. 

The SFR for IR detected galaxies was derived following 
Kennicutt (1998, see also Bell et al. 2005) normalized to a Kroupa 
(2001) IMF: 

SFR w+ „ = 1.15 x 10- 10 (L„ + 3.3 x L uv ) [M 0 yU‘] (1) 

where L uv is the luminosity at 280 nm rest-frame in erg s -1 . 

In order to estimate SFRs for galaxies not detected in the 

- The central wavelengths (widths) in nm of these filters are 359.84 (58.36), 
549.02 (85.79), 1237.59 (169.48) for the U, V, J bands, respectively. 
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redshift UV slope (/?) 

Figure 1. SFR ir versus redshift (left panel ) and IRX-/? relation (right panel ) for the sample of massive (M > 10 10 M©) galaxies at 1.0 < z < 1.5 detected by 
MIPS at 24 pm. The blue line is the corresponding SFR jr for the 24 pm MIPS detection limit (30 pjy at 5cr). The red points (23 galaxies) are IR-faint emitters, 
i.e., galaxies with an IR emission lower than the 5<x limit, selected to derive the IRX-/? relation. To avoid biasing the fits in the IRX-/3 plot, we randomly chose 
a subsample of faint IR emitters homogeneously covering the whole range of UV slopes /?. We fit the IRX-/? (log linear) relation for the faint IR emitters (red 
line) to obtain a /?-based dust attenuation for the population of galaxies with low levels of star formation. The green line is the IRX-/3 fit from Meurer et al. 
(1999), based on the analysis of regular star-forming galaxies and starbursts. Note that the dust attenuation correction at high /? values is much smaller for the 
IR-faint sample compared to IR-bright galaxies. Typical uncertainties are shown as error bars in the right corner of each panel. 


mid- or far-IR, we used UV luminosities at 280 nm rest-frame 
alone. These luminosities were converted into SFRs by applying 
Kennicutt (1998) equation normalized to a Kroupa (2001) IMF: 


SFRyv' = 0.98 x 10- 28 x L uv [M e yU 1 ] (2) 

The UV-based SFRs were corrected for dust attenuation fol¬ 
lowing a recipe based on the relation between the UV slope, /?, and 
the UV/IR ratio, IRX, know as the IRX-/? relation. The UV slope 
for each galaxy is calculated using a linear interpolation between 
150 and 280 nm in the templates fitting the SED of each galaxy 
(from which a photometric redshift and stellar mass estimate were 
obtained). The typical uncertainty in the /? values is ~ 20%. In or¬ 
der to convert these slopes to a dust attenuation, we used the com¬ 
parison between the observed UV and IR-based SFRs for galaxies 
detected by MIPS (Figure 1). Typical attenuation recipes based on 
the UV slope are derived for galaxies with high levels of star for¬ 
mation (Meurer et al. 1999; Takeuchi et al. 2012), which are dustier 
than relatively quiescent galaxies. Quiescent galaxies have system¬ 
atically lower ratio of total far-IR to UV luminosity than starburst 
galaxies (Kong et al. 2004). Given that in this paper we are inter¬ 
ested in galaxies with very low levels of star formation, and in order 
to avoid an over-correction for dust attenuation, we derived a IRX- 
/? relation for faint-IR emitters, i.e., those galaxies for which their 
IR detection is below the 5cr IR detection limit at that redshift. In 
comparison with bright-IR galaxies or the local starbursts, this sub¬ 
sample should present dust attenuation properties which are closer 
to the MQGs that we want to study. The procedure is outlined in 
Figure 1 , where we show SFR ir versus redshift, and the IRX-/? re¬ 
lation for the galaxies at 1.0 < z < 1.5 detected by MIPS. We high¬ 


light the galaxies which are faint-IR emitters, for which we derive 
the following IRX-/? relation: 

SFR 

IRX =-— = 8.09 + 3.02 x /? (3) 

SFR W 

We apply the Meurer et al. (1999) IRX-/? relation for/? values 
lower than the cross-point QS=-0.97) and Eq. 3 for higher /? values. 

2.2 Selection of quiescent galaxies 

To construct a complete and uncontaminated sample of MQGs at 
z=1.0- 1.5 we used two complementary methods: a UVJ diagram 
and sSFRs. We only consider galaxies with M > 10 10 M 0 , for which 
our survey is complete within the considered redshift range. 

2.2.1 Quiescence criterion based on the UV J colour-colour 
diagram 

There are 410 galaxies with M > 10 10 M Q at z=1.0 - 1.5 in the 
Rainbow catalogue of the GOODS-N field. The median stellar mass 
and sSFR of this sample are M ~ 10 104 M 0 and sSFR ~ 0.6 Gyr -1 . 
In Figure 2, we plot a UV J diagram including all these galaxies. 
Whitaker et al. (2011) defines the following as the region in the 
UV J diagram where quiescent galaxies are located: 

(U - V) > 1.3 
(V-J) < 1.6 

(U - V) > 0.875 x (V - J) + 0.6 (4) 

We find 87 galaxies with M > 10 1(> M 0 in the quiescent re¬ 
gion of the UVJ diagram. Out of these, 25% are detected by MIPS, 
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Figure 2. Plots showing the criteria to select the MQG sample. Left panel: UVJ diagram for galaxies at z= 1.0 - 1 .5 (grey dots). The final sample of MQGs (104 
objects) are marked with black empty circles. The continuous line delimits the region where quiescent galaxies should be located, according to Whitaker et al. 
(2011). The dashed line splits the previous region in two, with post-starburst galaxies located on the left. The f/W-selcctcd galaxies are coloured in green, 
and the sSFR-selected in red. The large blue circles mark IR detected galaxies (only for galaxies with sSFR < 0.2 Gyr 1 ). Right panel: sSFR versus stellar 
mass, colour coded as in the left panel. The dashed line represents our limit to consider a galaxy as quiescent (sSFR < 0.2 Gyr ! ). Uncertainties in the UVJ 
colours calculated by propagating the redshift uncertainties are negligible, given the superb quality of our SFlARDS-based photometric redshifts. Therefore, 
we assumed for our rest-frame UVJ absolute magnitudes the average uncertainty of the two observed biters which lie closer to the central wavelength probed 
by the UVJ rest-frame biters (shown as error bars in the right comer of the left panel). 


implying either some level of (obscured) star formation or nuclear 
activity. We eliminate them from this UVJ- selected sample of qui¬ 
escent galaxies. The sample of MQGs selected with the UVJ dia¬ 
gram and no IR detection is composed of 65 galaxies. The median 
and quartile stellar mass for this sample is log(M/M 0 ) = 10.7 Jq ®, 
they lie at z = 1.17}^, and have sSFR= 0.07^3 Gyr -1 . We will re¬ 
fer to the galaxies selected in this way as “UVJ-selected” galaxies 
and we will discuss their properties in Sect. 4.2. 


2.2.2 Quiescence criterion based on sSFR values 

We complement the selection based on the UVJ diagram with a cut 
in sSFR. We arbitrarily choose sSFR < 0.2 Gyr -1 as our limit for 
quiescence, given that by imposing this value we are able to virtu¬ 
ally recover all the sources identified as dead by the UVJ criterion 
(see right panel of Figure 2). 

We find 102 galaxies with M > 10 10 M 0 and sSFR < 0.2 Gyr -1 . 
All galaxies identified as quiescent in the UVJ diagram have 
sSFR < 0.2 Gyr -1 , except 2. These 2 galaxies have relatively 
young stellar populations (t < 0.8 Gyr) and short star formation 
timescales (r < 60 Myr), as we discuss in Sect. 3. Note that 
among the sSFR-based sample we have 22 IR emitters (21% of the 
sSFR-selected sample). These galaxies must have some (residual) 
star formation or nuclear activity, but their mass is large enough 
(log(M/M 0 ) > 10.4) to present sSFRs values comparable to dead 
galaxies undetected in the IR. The median and quartiles for the stel¬ 
lar mass distribution of this sample is log(M/M 0 ) = 10.8M q , 
they lie at z = 1.15} ^, and have sSFR= 0.07® }® Gyr -1 . With the 
sSFR criterion, we recover 8 galaxies located in the UVJ quiescent 
region which were discarded in the colour-colour selection because 
of their IR detection. However, they have low enough sSFR values 
(< 0.18 Gyr -1 ) to be selected in our MQG sample. We will refer to 
galaxies selected in this way as “sSFR-selected” galaxies, and we 
will discuss their properties in Section 4.2. Note that this is a com¬ 


plementary sample to the (/Vi-selected, i.e., they are galaxies with 
sSFR < 0.2 Gyr -1 but excluding the 65 previously selected with the 
UVJ colour criteria, which results in 39 additional galaxies. 

Combining the UVJ and sSFR criteria, we arrive to a final 
sample of MQGs with M > 10 10 M 0 composed of 104 sources. In 
Figure 3 we show the postage stamps made with HST ACS/WFC3 
for this final sample of MQGs. 


2.2.3 Statistical properties of the sample 

The combined UVJ- and sSFR-selected sample of MQGs was re¬ 
fined by carrying out a detailed stellar population synthesis analysis 
of each galaxy. We constructed the most detailed SED possible for 
each galaxy combining the ultra-deep SHARDS data with the G102 
and G141 spectro-photometric observations. To increase the S/N of 
the grism spectra, we binned them in order to have 10 nm per pixel. 
This corresponds to one and two resolution element for G102 and 
G141, respectively. This provides SEDs with up to 150 photomet¬ 
ric points at a photo-spectral resolution from 500 nm to 1700 nm. 
Galaxies without G141 or G102 spectra have at least 30 photo¬ 
metric data points. This unique photometric dataset encompasses, 
within the whole observed redshift range, significant spectral fea¬ 
tures related to the age of the galaxies, such as D4000 or Mgyy. 
Our final sample of MQGs at z = 1.0 - 1.5 consists of 104 galax¬ 
ies, 65 {/VT-selected plus 39 sSFR-selected. There are 54 (52%) 
galaxies for which spectroscopic redshifts are available, and the 
photometric redshift quality for them is characterized by a median 
Az/(1 + z) =0.0047. They are detected at 3cr in at least 13 SHARDS 
bands. Concerning the availability of grism data, 60% and 70% of 
the final sample have usable G102 and G141 spectra with at least 
S/N ~ 3 and median S/N ~ 10 per pixel (i.e., 10 nm). The rest have 
either severe contamination problems or are too faint for the grism 
observations. MQGs represent ~ 25% of the population of massive 
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Figure 3. 5"x5"postage stamps made with HST ACS/WFC3 for the 104 MQGs selected in this work (there are 2 sSFR-selected galaxies missing because 
they are outside the WFC3 region). Galaxies in the upper panel are f/VJ-selected, while galaxies in the lower panel are sSFR-selected. In each panel, they are 
ranked (from left to right, top to bottom) by increasing mass-weighted age t m (see Sect. 4). Postage stamps for the 3 sub-samples of galaxies defined in the 
text based on t\j are framed in blue (mature galaxies), green (intermediate), and red (senior galaxies). 


galaxies within the same redshift interval. They have median values 
log(M/M s )= 10.7“;° M 0 , z = and sSFR= 0.10°;^ Gyr 1 . 

The fraction of quiescent galaxies in a mass-selected sample 
according to our analysis is in good agreement with the 28% re¬ 
ported by Muzzin et al. (2013) for a purely UVJ selected sample 
of quiescent galaxies with log(M/M 0 ) > 9.5 at z = 1.0 — 1.5. The 
fraction of MQGs that we find is larger than the 13% of quiescent 
galaxies with log(M/M Q ) > 9.6 at z = 1.1 - 1.5 selected on the ba¬ 


sis of their NUV, r + and J colours found by Ilbert et al. (2013). The 
fraction of very massive (log(M/M 0 ) > 10.85 M 0 ) quiescent galax¬ 
ies (selected on the basis of IR colours and sSFR sed < 0.01 Gyr -1 ) 
derived in Fontana et al. (2009) at z ~ 1.2 is ~ 40%, larger than 
the ~ 25% that we obtain with the same mass cut. We will further 
discuss our results about the properties of our sample of MQGs and 
compare them with the literature in Sect. 4.3. 
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Figure 4. Example of SED-fitting for one of our MQGs, SHARDSJ123727.9+622034.7, a galaxy at z =1.148. In the left panel we show the SED-fitting to the 
whole dataset gathered for our work (194 data points). The red line is the best-fitting model. The white circles represent the SHARDS spectro-photometric data 
and the diamonds are broad-band data. The G102 and G141 spectra are plotted as black and dark grey lines and zoomed in the central and bottom right panels. 
We depict 3cr errors in all plots (and the average S/N in the grism spectrum plots). The vertical dotted lines show the location of typical emission lines (Her, H/? 
and [fill [,13727). We also show a 5"x5"postage stamp made with HST ACS/WFC3 data. In the upper right panel, we show a zoom in the SHARDS region. 
The coloured areas represent the bands used to determine the Mgyv and D4000 indices, whose values are given in the legend. For this galaxy, only one stellar 
population model was compatible with the data after applying the method described in Sect. 3.2. The best-fitting parameters are shown in the legend, including 
metallicity, star formation timescale, age, mass-weighted age, dust attenuation, stellar mass, predicted 24 //m flux (see Sect. 3.3) and statistical significance of 
the solution. 


3 METHODOLOGY TO DETERMINE THE STAR 

FORMATION HISTORIES OF MASSIVE QUIESCENT 
GALAXIES 

3.1 SED-fitting 

Our main goal in this paper is to study in detail the SFHs of MQGs 
at 1.0 < z < 1.5. For this purpose, we fit the observed photome¬ 
try to stellar population synthesis models using the synthesizer fit¬ 
ting code described in Perez-Gonzalez et al. (2003, 2008). We use 
models with one burst of star formation characterized by a delayed 
exponentially declining SFH: 

SFR(t) oc f/r 2 x e~" T (5) 

The most common SFH parameterization used in the litera¬ 
ture is an exponentially declining function. However, we have cho¬ 
sen a more realistic parameterization with an initial increase in the 
star formation activity followed by an exponential decay, avoid¬ 


ing the nonphysical infinite derivative at time equals zero obtained 
in the pure exponential. Our parameterization is also closer to 
the SFHs predicted by galaxy evolution models (e.g. Pacifici et al. 
2013; Barro et al. 2014). 

We compare our SEDs with the BC03 stellar population li¬ 
brary, assuming a Kroupa (2001) IMF and the Calzetti et al. (2000) 
attenuation law. We allow 3 different metallicity values, Z/ZO = 
0.4, 1.0 and 2.5, which correspond to values of sub-solar, solar 
and super-solar metallicity. Although it is commonly assumed in 
the literature that galaxies have solar metallicity, it is well known 
that galaxies may have different metal contents. For example, the 
mass-metallicity relation (Tremonti et al. 2004) points out that, at 
low redshift, more massive galaxies are more metal-rich than less 
massive ones. The metallicity also affects the shape of the SED, 
making galaxies look redder as we move to higher metallicities. 
Although the metallicities may take on more values than assumed, 
using only the 3 discrete values around solar metallicity given by 
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Parameter 

range/values 

units 

step 

age (f 0 ) 

0.04-6.3 

Gyr 

0.1 dex 

timescale (r) 

3 - 10000 

Myr 

0.1 dex 

dust attenuation (A v ) 

0-1.5 

mag 

0.1 mag 

metallicity (Z) 

0.4, 1.0, 2.5 

z/z © 

discrete 


Table 1. Free parameters (age, star formation timescale, dust attenuation 
and metallicity) and their allowed ranges used in the SED-fitting procedure. 

the BC03 models is a sufficient approximation to see the effect of 
this parameter in our results. 

Our stellar population synthesis code, synthesizer, performs a 
X 1 minimization and returns the model which best fits the data. The 
age (time since the formation of the stellar population that domi¬ 
nates the SED of the galaxy, to), star formation timescale (r), dust 
attenuation (A v ) and metallicity (Z) are set as free parameters. The 
parameter space allowed in the fitting procedure is given in Table 1. 
The stellar mass (M) of each galaxy is derived as the normalization 
of the best-fitting template to the median observed flux. 

In the following sections, apart from the stellar population 
parameters mentioned above, we will also discuss mass-weighted 
ages (t M ). The t M is defined as: 

r°SFR(t) x t df C° SFR(t) x t dt 

t M = - = -- (6) 

J^°SFR(t)df M 

where to is the best-fitting age, and M is the best-fitting mass. While 
to corresponds to the beginning of the star-formation period that 
dominates the SFH of the galaxy, t M is a better approximation to 
the average age of the stellar population and is a more useful pa¬ 
rameter to understand the evolution of galaxies taking into account 
the duration of the star formation processes. 

An example of the SED-fitting is presented in Figure 4, where 
we show the fit to the whole wavelength interval, as well as zooms 
in the spectral regions covered by SHARDS, and the G102 and 
G141 grisms. Our fitting code also includes an algorithm to es¬ 
timate uncertainties in the derived parameters and to study (and 
break when possible) the degeneracies typically present in SED- 
fitting studies. This algorithm involves a Monte Carlo technique, 
and the usage of direct measurements of spectral indices (especially 
D4000 and Mg uv ) with the SHARDS and WFC3/grism spectro- 
photometric data, as well as the analysis of the mid- and far-IR 
fluxes and upper limits. We describe this algorithm in detail in the 
following subsection. 

3.2 Estimating uncertainties and analysing degeneracies 
with a Monte Carlo algorithm 

We used a Monte Carlo approach to estimate uncertainties in the 
stellar population properties and to take into account the possi¬ 
ble degeneracies of the SED-fitting technique. For each galaxy, 
we constructed 1000 modified SEDs by allowing the photomet¬ 
ric data points to randomly vary following a Gaussian distribution, 
with a width given by the photometric errors. We performed the 
SED-fitting to the modified photometric data and obtained 1000 
different solutions with their corresponding set of parameters for 
every galaxy in our sample. In the SED-fitting procedure, the stel¬ 
lar population modeling code looks for best-fitting ages, timescales, 
metallicities, and dust attenuations in a grid of discrete values. 

Once we have 1000 SED-fitting solutions for a given galaxy, 


we look for clusters of solutions in the r-t 0 parameter space. In 
order to account for the discrete distribution of the fitting parame¬ 
ters probed by the minimization algorithm, we introduce Gaussian 
noise to the output parameters of each galaxy using a width equal 
to the step used for each fitted property (see Table 1 ). We identify 
clusters in the r-to plane with a k-means method and a minimal 
separation of 0.2 dex between different solutions (i.e., the differ¬ 
ence between the median cluster properties must be at least 0.2 dex 
in age and r). Solutions which provide similar results are grouped 
as a single solution identified by a median value and a scatter in the 
multi-dimensional t 0 -r-A v -Z space. Using the full set of solutions 
for a given cluster, we calculate the values enclosing 68% of the 
data around the median. These are assumed to be the uncertainties 
of our estimations for each cluster of solutions. The typical rela¬ 
tive uncertainties of the parameters in our analysis are Ato=12%, 
Ar=16%, AA V =0.06 mag and AM=0.05 dex. The metallicity un¬ 
certainties cannot be determined because the allowed Z values are 
discrete. 

Each cluster is assigned a statistical significance given by the 
fraction of solutions belonging to that cluster. The clusters repre¬ 
sent the significant solutions of a galaxy taking into account the 
degeneracies in the determination of the stellar properties from the 
SED-fitting. Although we use the t 0 -r plane to look for different 
solutions, we tested whether looking for clusters in the whole to- 
r-A v -Z multi-dimensional space changed the results. Given that all 
the parameters are highly correlated, we found no difference be¬ 
tween both approaches, i.e., a cluster analysis in one plane was able 
to robustly recover clusters in the 4-dimensional space. 

In Figure 5, we show the SED-fitting result for a galaxy for 
which we find three clusters of solutions. The most significant (74% 
of solutions belong to this cluster) is a relatively young stellar pop¬ 
ulation (t M ~ 1.4 Gyr) with star formation timescale r ~ 160 Myr, 
moderate dust attenuation (A v = 0.3 mag), and super-solar metal¬ 
licity. Another solution consistent with the data is an older pop¬ 
ulation ( t M ~ 2.6 Gyr) with a longer star formation timescale 
(r ~ 310 Myr), higher dust attenuation (A v = 0.6 mag), and sub¬ 
solar metallicity. And finally, another possible solution is charac¬ 
terized by an intermediate age population {tM ~ 1.9 Gyr) with very 
short star formation timescale (r ~ 24 Myr), higher dust attenua¬ 
tion (A v = 0.7 mag), and also sub-solar metallicity. We remark that 
all 3 solutions fit the data equally well, although we must warn the 
reader that most of the spectro-photometric data points are bluer 
than 1 ftm rest-frame, so the fits are biased towards the bluer part 
of the SED. In Figure 6, we show the procedure used to identify 
the three clusters of solutions. We plot in the to-r space the solu¬ 
tions obtained for the 1000 Monte Carlo simulations. The 3 clusters 
mentioned above can be identified by colours, and we also depict 
the median values of each solution. 

We investigate how our results are affected by the degenera¬ 
cies by analysing the clusters of solutions. We obtain 176 different 
solutions for the 104 galaxies in our sample. This means that each 
galaxy has, on average, 1.7 solutions. There are 48 galaxies (~ 46% 
of the sample) for which we find only one cluster of solutions, i.e., 
only one set of properties fits the data, after taking into account 
observational uncertainties and the degeneracies linked to them. 
There are 41 galaxies (39%) with 2 clusters of solutions and less 
than 15% of the sample has 3 or more clusters (with the maximum 
number of clusters being 4, which happens only for two galaxies). 
The fact that almost half of the sample presents no degenerate so¬ 
lutions reveals the power of combining SHARDS and WFC3 data 
to constrain the properties of MQGs at high-z. When using broad- 
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Figure 5. Example of our SED-fitting method, including the analysis of degeneracies, for SHARDSJ 123620.3+6200844.3, a galaxy at z sp = 1.016. Colours 
and symbols as in Figure 4. For this galaxy, three clusters of solutions were obtained after applying the method described in Sect. 3.2. The main parameters of 
each solution are shown in the legend. 


band data alone, the number of galaxies with degenerate solutions 
increases up to 76%. 

We derived the maximal difference of the galaxy parameters 
obtained from the different clusters of solutions for each galaxy. 
The mean relative differences are (5to ~ 40%, St M ~ 30%, St ~ 80%, 
(5A v =0.1 mag, <5M ~ 0.1 dex. The differences between solutions 
for the stellar masses and the dust attenuations are of the order of 
the typical uncertainties of each cluster of solutions, meaning that 
these parameters are not strongly affected by the degeneracies. In 
fact, only two galaxies with degenerate solutions present (SA V val¬ 
ues larger than 0.5 mag. The degeneracies in to are smaller than 
50% in relative values (and even smaller (30%) for the t M ). The 
largest differences are found for the star formation timescale, r. 
Again, the metallicity uncertainties cannot be determined because 
the allowed Z values are discrete. The metallicity values are unique 
(i.e., SZ =0) for 55% of the galaxies with more than one cluster of 
solutions. 


3.3 Using spectral indices and IR detections to break 
degeneracies 

The SED-fitting technique described in the previous section does 
not make use of the full power of the dataset gathered for this 


work. Including all photometric data points from the UV to the 
mid-IR gives us the global shape of the stellar emission. But this 
global shape may easily wash out the higher spectral informa¬ 
tion given by the spectro-photometric data from SHARDS and the 
grism observations. Thus, the stellar population properties can be 
even better constrained when taking advantage of the ultra-deep 
spectro-photometric SHARDS and grism data. The high resolution 
(R ~ 50) photometry from SHARDS and the grisms allows us to 
measure spectral indices such as Mguv and D4000, which are well 
correlated to stellar population properties (Kauffmann et al. 2003; 
Daddi et al. 2005). 

We measure the Mguv index presented in Daddi et al. (2005) 
as a strong age-dependent feature in the UV, detectable in relatively 
low resolution spectra in the region 260-290 nm. This feature is an 
absorption band formed by a combination of several strong Mg and 
Fe lines. The Mguv index is defined as: 


(7) 


where the integration ranges are defined in nm. Note that to mea¬ 
sure the Mguv index we need 10 nm windows at rest-frame, which 
at redshift z > 1 translates to 20 nm or more (observed-frame). 
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Figure 6. Example of the procedure used to break degeneracies among best-fitting solutions using measurements of absorption indices for 
SHARDSJ123620.3+6200844.3. Upper panel : we show a zoom of the SED, centered in the SHARDS region, for the same galaxy depicted in Figure 5. 
The shaded bands mark the regions used to measure the Mguv and D4000 indices. The coloured circles represent the data points used to measure the indices 
from the SHARDS photometry. Lower left panel'. Distribution of the 1000 best-fitting solutions in the age/timescale plane. Note that this is the direct age from 
the SED-fitting (to) and not l\r. The coloured dots correspond to each one of the identified clusters, with the same colour as in the SED plot. For each one 
of the three clusters, the median values are represented by large green circles and the black contours enclose 68% of the solutions. Lower right panel', we 
show the evolutionary tracks in the Mgr/y-D4000 plane for the three best-fitting models found for SHARDSJ123620.3+6200844.3. The empty circles show 
the expected values of the indices at different ages, to (from smaller to larger symbols: 0.5, 1.0, 2.0, 3.0 and 5.0 Gyr). The large coloured diamonds show the 
location of the indices at the age of the best-fitting solution. The black circle represents the indices measured from the SHARDS photometry. In this case, 
unlike in the following plots where we will use tM, to is the appropriate parameter to separate clusters of solutions and to identify the evolutionary tracks 
consistent with the indices, as to is the actual best-fitting parameter of the stellar population models. The indices predicted by the track of the first solution are 
incompatible with the values measured directly from photometry and therefore, this solution is discarded. 


This means that the SHARDS filters, with a typical width around 
16 nm, present a sufficient spectral resolution to measure the Mguv 
absorption index. 

We also measured the D4000 index, introduced by Bruzual 
(1983), which gives an estimate of the strength of the 4000 A break. 
This index is defined as: 



where the integration ranges are again defined in nm. 

We prefer to use this definition instead of the narrow index 
D„4000 (Balogh et al. 1999) to reduce the uncertainty measure¬ 
ment of the index using photometric data alone, thanks to the 
broader index bands. See Hernan-Caballero et al. (2013) for more 
details about measuring D4000 and D„4000 with SHARDS data. 


To measure the spectral indices based on our spectro- 
photometric data, we first select the data points for each galaxy that 
fall into the bands used in the index definition. These bands depend 
on the redshift of the galaxy, but we note that for SHARDS we also 
have to consider the position of the galaxy in the FOV, given that the 
central wavelength of the pass-band seen for each galaxy depends 
on the position in the OSIRIS focal plane (see Perez-Gonzalez et al. 
2013). We directly measured the ratios between the corresponding 
fluxes to get a first estimation of the indices. We refer to these in¬ 
dices measurements as Mguv* and D4000*. To be able to compare 
these values with the typical index definitions, we correct them by 
using the ratio between the index measured in the best-fitting model 
of each galaxy at the central wavelength of each filter convolved to 
the SHARDS resolution and the standard index measured directly 
in the models. The typical value of this correction is small, the av¬ 
erage is 1.11 for Mg uv , and 0.99 for D4000. We also checked that 
the correction is rather insensitive to the use of only one stellar pop- 
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ulation model to measure the correction, i.e., using a larger set of 
models with different timescales and ages has a very small effect (< 
10%). We remark that the indices calculated in this way are com¬ 
pletely independent of the SED-fitting procedure and are based on 
the observed photometric data alone. 

In Figure 6, we plot a zoom of the SED in the spectral range 
covered by SHARDS for the galaxy presented in Figure 5. We 
also plot the corrected Mg uv and D4000 indices for that partic¬ 
ular galaxy, as well as the tracks expected for the evolution of these 
indices according to the models of the 3 best-fitting (degenerate) 
solutions obtained for that galaxy (see Sect. 3.2). The indices mea¬ 
surements allow us to eliminate solutions which are incompatible 
with the spectro-photometric data. For this particular galaxy, the 
measured indices are more compatible with the second and third 
solutions, and we discard the most significant solution (the one with 
to ~ 1.7 Gyr). 

In addition to the use of spectral indices, we evaluated the 
robustness of the different solutions obtained for each galaxy by 
making use of the (un-)detection in the mid- and far-IR and energy 
balance arguments. For each best-fitting solution, given the dust at¬ 
tenuation and mass, we derived the expected flux at an observed 
wavelength of 24 yum (F pr (24), for the galaxy shown in Figure 4, 
the values are written in the legend). The procedure starts calcu¬ 
lating the luminosity absorbed by dust in the UV/optical, which 
has to be re-emitted in the IR. We assumed that the stellar emis¬ 
sion absorbed by dust must be equal to the IR luminosity integrated 
from 8 to 1000 /jm, L/ R . We then used the Rujopakarn et al. (2013) 
relation to transform the h /R into a 24 [i m observed flux (a rela¬ 
tion that depends on redshift). We note that the Rujopakarn et al. 
(2013) templates are based on star-forming galaxies, which may 
not be comparable to the quiescent galaxies from this sample. 
However, we consider it a good approximation, since the contri¬ 
bution to the dust heating of the older stellar population is larger 
at wavelengths > 250 [i m (see Bendo et al. 2012) and does not 
significantly affect the 24 fjm emission.’ For galaxies undetected 
in the IR, we were able to eliminate solutions which present ex¬ 
pected 24 /jm fluxes larger than the detection limit (~ 50 /jJy at 
50% completeness, as estimated in Perez-Gonzalez et al. 2005 to 
build the IR luminosity function at z = 0 — 3). For example, for 
the SHARDSJ123620.3+6200844.3, the third solution predicts a 
24 [i m flux of ~ 55 /rJy. As this galaxy is not detected in the IR, 
we rejected the third solution and chose as the best solution the 
second one ( t M ~ 2.6 Gyr, r =300 Myr). The second solution has 
F pr (24) — 51, which is still larger than the 50% completeness limit, 
but only by 1%. With this method, we were able to reject only ~ 4% 
of the degenerate solutions, but it was a strong argument to limit the 
dust attenuation range probed by our stellar population analysis to 
values within 0 < A v < 1.5 mag. Larger values would imply IR 
detections for galaxies as massive as ours. 

We are able to break the degeneracies making use of the spec¬ 
tral indices for ~ 32% of the galaxies with more than one cluster 
of solutions. For the remaining objects, the indices uncertainties 
are too large (due to photometric errors) or the indices predicted 
by the tracks of the different solutions are consistent with the mea¬ 
sured values. The energy balance argument helps breaking the de- 


3 Since for the same L/«, the predicted 24 jjm flux for a quiescent galaxy 
should be lower than for a star-forming galaxy, our predicted 24 fjm fluxes 
calculated with the method described above may be overstated. This reduces 
the significance of our rejection of non IR detections, although we only do 
so for 4% of the degenerate solutions. 


generacies in 4% of the cases. When all the solutions of a galaxy 
are compatible with the measured indices and do not violate the 
energy balance argument, we choose the most significant one, i.e., 
the most populated cluster (~ 32% of the cases). For 22% of the 
galaxies, the solutions were very similar in t - A, - Z, with only 
significant differences in r. We also discarded 10% of the solutions 
for being unrealistic (r < 10 Myr and log(M/M 0 ) > 11.0, which 
would imply SFR ~ 10000 M 0 yr~* or larger). 

In summary, out of the complete sample of 104 galaxies, 46% 
of them had only one possible solution (i.e., no degeneracies). Out 
of the 54% of galaxies with degenerate solutions, we were able 
to break the degeneracies either by measuring indices or by us¬ 
ing the Mid-IR/Far-IR data for 20% of the sample. For the remain¬ 
ing ~ 34%, we used the most significant solution. In the following 
sections, we will only consider one solution for each galaxy and we 
will refer to them as primary solutions. In only 12% of the cases the 
primary solutions are not the most significant ones. The properties 
of the primary solutions for each galaxy are given in Table 2. 


4 ANALYSIS OF THE STAR FORMATION HISTORIES 
OF MASSIVE QUIESCENT GALAXIES 

In this section, we analyse the stellar populations properties of 
MQGs at 1.0 < z < 1.5. For the discussion, we only consider 
the primary solutions identified with the methodology described in 
Section 3. 

4.1 Statistical properties of the stellar populations of MQGs 

In Figure 7, we show the location of our galaxies in the mass- 
weighted age versus mass (t^-M) plane, colour-coded using 3 
bins in star formation timescale. The r bins have been chosen at 
the values were the t M -T relation changes trend, according to the 
right panel of the same figure, where we plot the relationship be¬ 
tween t M and r, using a colour-code to distinguish 3 bins of mass. 
The mass bins limits are approximately the median and 3 rd quar- 
tile values. The population of MQGs at z = 1.0 - 1.5 are domi¬ 
nated by “new arrivals”, i.e., galaxies with relatively young stellar 
populations: t M < 2 Gyr. These galaxies younger than 2 Gyr ac¬ 
count for 85% of the sample. We also identify a tail of old galax¬ 
ies ( t M > 2 Gyr) summing up ~ 15% of the total population. 
Hereafter we divide the MQGs in three sub-samples depending on 
their t M values: mature ( t M <1.0 Gyr, 38 galaxies), intermediate 
(t M =l- 0-2.0 Gyr, 50 galaxies), and senior ( tM > 2.0 Gyr, 16 galax¬ 
ies). The statistical properties of the galaxies divided in age and 
mass bins are listed in Table 3 and 4, respectively. 

The population of mature galaxies presents an average mass- 
weighted age (t M )= 0.8q7 Gyr and relatively short timescales 
(t)=602q° Myr. The typical mass is log((M)/M o )=10.4|°3. On 
the other hand, the senior population presents average values of 
(tM)=2.62 2 Gyr, longer star formation timescales, (t)=400jqq Myr, 
and larger masses log((M)/M Q )=10.7[g°. The intermediate popula¬ 
tion has transitional parameters: {t M )= 1.4] 3 Gyr, (r)=2003°° Myr 
and log«M>/Me)= 10.5}“. 

Note that the fact that the mature population has short r values 
is a direct consequence of our sample selection. As we are select¬ 
ing quiescent galaxies (galaxies with low levels of sSFR), galaxies 
with longer r (> 100 Myr) and young ages (< 1 Gyr) would have 
too high sSFR values to enter in our quiescent selection criteria. In 
principle, the senior galaxies might present short or long star for¬ 
mation timescales, i.e., we do not have any selection bias against 
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ID 

z 

V-J 

U-V 

sSFR ir+uv 

log(M/M 0 ) 

to 

tM 

T 

A v 

Z/Zo 

SFRsed 

Sign. 





[Gyr 1 ] 


[Gyr] 

[Gyr] 

[Myr] 

[mag] 


[M 0 yr->] 

% 

SHARDS 123737.94+621309.0 

1.2410 (s) 

1.34 ± 0.09 

1.95 ± 0.09 

0.10 ±0.06 

10.78 ± 0.05 

2 l 2 - 4 

Z -M.9 

i-7?:2 

199 223 

0.00 ± 0.02 

2 5 2-5 

Z ' J 2.5 

0.08 

100 

SHARDS 123657.46+621451.2 

1.2534 (s) 

1.20 ± 0.08 

2.00 ± 0.07 

0.02± 0.01 

11.09 ±0.05 

2 ' 4 22 

L8 U 

3 1 6 354 
J1U 280 

0.66 ± 0.07 

0 4 0 - 4 
T).4 

1.40 

100 

SHARDS 123723.91 +621520.7 

1.39 (p) 

1.0 ±0.2 

1.7 ±0.2 

0.1 ±0.1 

10.0 ±0.1 

°- 8 n? 

°-C 

14 35 

0.2 ± 0.3 

2 5 2-5 

< nr 3 

100 


Table 2. Galaxy properties for the sample of 104 galaxies here presented: ID, redshift (z; s for spectroscopic, p for photometric), UVJ colours, sSFR used in 
the sample selection (see Sect. 2.2), mass (M), best-fitting age (to), mass-weighted age (Im), star formation timescale (r), dust attenuation (A v ), metallicity 
(Z), SFR from SED-fitting (SFRsed) and significance of the primary solution. The full version of the table is available online in the supplementary files. 
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Figure 7. Mass-weighted ages versus stellar mass (left pane!) and star formation timescale (right panel). Galaxies are colour coded by star formation timescale 
and mass, respectively. Only the primary solutions (those selected with the different methods described in Sect.s 3.2 and 3.3) are plotted. Circles represent 
galaxies with only one solution, while triangles represent galaxies with degenerate solutions. The light red areas represent the region between the l If and 3 rd 
quartiles for each parameter. Diamonds represent the median values for the whole sample (yellow) and for the subsamples (colour coded according to the 
legend of each panel). The coloured error bars show the region encompassing 68% of the data for each subsample. Typical individual uncertainties are plotted 
in the lower right corner of each plot. 


any of those types of galaxies. The bias in our results due to our 
sample selection are further discussed in Sect. 4.3, together with 
the average SFHs of galaxies divided in mass and age bins. 

With respect to the dust attenuation, the senior galaxies are 
less dusty, (A v )=0.4jJ j mag, than the mature ones, (A v >=0.8 ( l ) , j mag. 
The attenuation values are in agreement with those found by other 
authors (e.g., Belli et al. 2015). The galaxies with the highest dust 
attenuations are among the youngest (~ 90% of the galaxies with 
A v > 0.8 mag are younger than 2 Gyr), suggesting that they could 
be recently quenched galaxies. 

With respect to the metallicities, 41% of the galaxies are best- 
fitted with solar metallicity, 36% with super-solar metallicity and 
23% with sub-solar metallicity. No clear metallicity trends are 
found with respect to mass or age. 

At this point, we want to mention the difference between old 
and quiescence when using r-models. In single stellar population 
(SSP) models, all the star formation takes place at the same time 
(t 0 ), and then the population passively evolves. Therefore, an old 
galaxy is always quiescent and, in fact, a young galaxy would also 
be literally quiescent (although it may be not selected with our cri¬ 
teria). When using r models, the duration of the star formation 
depends on r. Galaxies with long star formation timescales may 
have been formed more than 1-2 Gyr ago but could still be form¬ 
ing stars. What actually indicates the quiescence of a galaxy (and 
consequently the star formation activity) is the sSFR or the t M -r 
relation. 


4 


3 


L? 2 


1 


0 

10 100 1000 
t [Myr] 

Figure 8. Mass-weighted ages versus star formation timescale. Galaxies are 
colour coded by their sSFRsed- Circles represent /7V 7 ./-selected galaxies, 
while the diamonds are sSFR-selected. The sizes of the symbols represent 
the percentage of mass assembled in the last Gyr. The grey areas represent 
the region between the 1” and 3 rd quartiles. Typical uncertainties are plot¬ 
ted in the lower right corner. 
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In Figure 8, we show again the t M versus r plot, but now the 
galaxies are colour-coded according to their sSFRsed (averaged 
over the last 100 Myr from the best-fitting models, i.e., these sSFRs 
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Parameter 

Sample 

Number 

Qi 

Median 

Q3 

Z 

Mature 

N=38 

1.12 

1.23 

1.33 


Interm. 

N=50 

1.02 

1.17 

1.24 


Senior 

N=16 

1.02 

1.06 

1.19 

to [Gyr] 

Mature 


0.8 

0.9 

1.0 


Interm. 


1.3 

1.8 

2.2 


Senior 


3.0 

3.2 

4.2 

r [Myr] 

Mature 


20 

60 

100 


Interm. 


30 

200 

300 


Senior 


300 

400 

500 

A v [mag] 

Mature 


0.5 

0.8 

1.1 


Interm. 


0.1 

0.5 

0.8 


Senior 


0.1 

0.4 

0.6 

Z/Z G 

Mature 


1.0 

1.0 

2.5 


Interm. 


0.4 

1.0 

2.5 


Senior 


0.4 

1.0 

1.0 

l°g(M/M 0 ) 

Mature 


10.3 

10.4 

10.7 


Interm. 


10.3 

10.5 

10.8 


Senior 


10.6 

10.7 

11.0 

hi [Gyr] 

Mature 


0.7 

0.8 

0.9 


Interm. 


1.2 

1.4 

1.7 


Senior 


2.2 

2.6 

3.4 

sSFRsed [Gyr" 1 ] 

Mature 


< 10 -4 

< 10 -4 

10 -2 


Interm. 


< 10 -4 

10 -3 

10 -2 


Senior 


10 -3 

10 -3 

10 -2 

t 9 [Gyr] 

Mature 


0.2 

0.5 

0.7 


Interm. 


0.7 

0.9 

1.1 


Senior 


1.5 

1.9 

2.6 


Table 3. Galaxy properties of our sample in three age bins: mature ( t m <1.0 
Gyr), intermediate (1.0 < tM < 2.0 Gyr) and senior (1m > 2.0 Gyr). We 
show the I s ', median and 3 rrf quartiles for redshift (z), the best-fitting ages 
(to), star formation timescales (r), dust attenuations (A v ), metallicity (Z), 
masses(M), mass-weighted ages Vm), sSFRsed and time since quenching 
(tq, explained in Sect. 4.2) 


are different from those used in the selection, which were based 
on the observed UV luminosity and dust attenuation estimates - 
see Section 2.2). The galaxies with similar sSFRsed values are dis¬ 
tributed along diagonal lines of constant t M lr values. The senior 
population of galaxies with long r values show higher sSFR val¬ 
ues than mature galaxies with very short r. The symbol size repre¬ 
sents the percentage of mass assembled in the last Gyr. The mature 
galaxies (? M <1 Gyr) assembled most of their mass during that time 
(by definition of mature objects). There is a population of inter¬ 
mediate galaxies with r >100 Myr that has assembled between 
10-50% during the last Gyr, while the senior population has assem¬ 
bled less than 10% of its mass. We also differentiate in Figure 8 
between l/W-selected and sSFR-selected galaxies. sSFR-selected 
objects, i.e., galaxies outside the quiescent (/VT-region (or inside 
the (7 W-quiescent region but detected in the IR, see Sect. 2.2), are 
preferentially located in the lower right corner of the t M -T plane. 
This indicates that the sSFR-selected galaxies have larger t^/r val¬ 
ues and are less quiescent than the pure UVJ galaxies. In Table 5 
we show the median properties of each sub-sample and we discuss 
the distribution of the derived galaxy properties in the i/W-plot in 
Sect. 4.2 and Figure 11. 

Our result about a duality (mature versus senior systems) in 
the population of MQGs at 1.0 < z < 1.5 and their differences 


Parameter 

Sample 

Number 

Qi 

Median 

Q3 

z 

Low M 

N=45 

1.13 

1.24 

1.33 


Interm. M 

N=36 

1.02 

1.14 

1.23 


High M 

N=23 

1.01 

1.12 

1.24 

to [Gyr] 

Low M 


0.9 

1.1 

1.6 


Interm. M 


1.0 

2.0 

2.3 


High M 


1.0 

2.0 

3.0 

t [Myr] 

Low M 


20 

80 

200 


Interm. M 


30 

200 

400 


High M 


30 

250 

300 

A v [mag] 

Low Mass 


0.2 

0.6 

1.0 


Interm. M 


0.3 

0.6 

0.9 


High M 


0.3 

0.5 

0.8 

z/z Q 

Low M 


1.0 

1.0 

2.5 


Interm. M 


1.0 

1.0 

2.5 


High M 


1.0 

1.0 

1.0 

log(M/M s ) 

Low M 


10.1 

10.3 

10.4 


Interm. M 


10.6 

10.7 

10.8 


High M 


10.8 

10.9 

11.0 

hi [Gyr] 

Low M 


0.8 

1.1 

1.3 


Interm. M 


1.0 

1.4 

1.8 


High M 


1.0 

1.6 

2.4 

sSFRsed [Gyr -1 ] 

Low M 


< 10 -4 

10 -4 

10" 2 


Interm. M 


< 10 -4 

10 -3 

10 -2 


High M 


< 10 -4 

10 -3 

10 -2 

t 9 [Gyr] 

Low M 


0.4 

0.6 

0.9 


Interm. M 


0.7 

0.9 

1.2 


High M 


0.8 

1.0 

1.7 

Table 4. Galaxy 

properties 

of our 

sample 

in 3 

mass 


bins: log(M/M o )=[10.0,10.5], log(M/M G )=[10.5,10.8] and 

log(M/M 0 )=[l().8,11.4], The parameters are the same as those in Ta¬ 
bles. 


not only in ages but also in timescales and dust attenuation was 
tested. The test was planned to check that the population of old 
(>2 Gyr) MQGs at z = 1.0 - 1.5, which is just a 15% and a tail in 
the age distribution, is real and its existence is not an effect of the 
degeneracies (in the T-t M plane). For this purpose, we constructed 
spectro-photometric stacks for the mature (t M <1.0 Gyr), interme¬ 
diate {t M - 1.0-2.0 Gyr) and senior (t M > 2.0 Gyr) populations. 

The broad-band colours of the three sub-samples of galaxies 
are very similar (they were selected with the same method based on 
colours) in the whole wavelength range. However, at the resolution 
achieved by the spectro-photometric data from SHARDS and the 
WFC3 grisms, especially in the Mguv and D4000 spectral regions, 
the sub-samples and the average fitting models show measurable 
differences. In Figure 9, we show the stellar population models 
with the average properties of the mature and senior populations, 
concentrating in the spectral region covered by SHARDS. There 
are significant differences in the relative fluxes for the two mod¬ 
els. In the D4000 region there is an excess of flux for the mature 
population model and the D4000 break is stronger for the senior 
population model. This translates to a ~ 20% difference in the rel¬ 
ative fluxes of the two models in the blue band of the D4000 index. 
In the Mguv region the relative flux of the mature model is even 
higher, up to 70%, revealing the signature of a younger population, 
which is less significant in the senior model. We should be able to 
distinguish between the two solutions by measuring the differences 
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Figure 9. Differences in the average stellar population models for the ma¬ 
ture and senior galaxy sub-samples (see text for details). Upper panel: stel¬ 
lar population models with the average properties of the mature and se¬ 
nior populations (as shown in the legend) in the spectral region covered by 
SHARDS. The models have been normalized to the fluxes in the 230-450 
nm range. We show both the original resolution of the models (thin lines) 
and the models convolved to the SHARDS resolution (R=50, thick lines). 
Lower panel: ratio of the two models shown above at the original model 
resolution (thin grey line) and at R=50 (thick black line). We remark the 
significant (and measurable) differences between the two models at the res¬ 
olution of our spectro-photometric data, especially for the the Mgpy and 
D4000 indices (marked in both plots as coloured areas). 

in the indices, and this is possible in the SHARDS and grism data 
(and not with broad-band observations). 

In Figure 10, we present the stacks built for the mature and 
senior galaxy populations together with the models characterized 
by their average properties. We also show the indices measured in 
the stacks of the mature, intermediate and senior galaxy popula¬ 
tions, and the tracks predicted by the model characterized by the 
average timescales and dust attenuations quoted in Table 3 for each 
sub-population. The indices measured directly from the stacks are 
consistent with those predicted by the tracks of the models with the 
average properties of each sub-population. We demonstrate that the 
spectro-photometric data directly show the differences in the stellar 
population properties for the 3 sub-samples. Note that the origin 
of the average stellar population properties is the full-SED spectral 
fits, and here we are comparing those models with finer resolution 
data in a limited wavelength range. The two analyses are not com¬ 
pletely independent, but the information encoded in the data at the 
resolution which makes index measurements possible could easily 
be erased or degraded by fitting the whole SED from the UV to 
the mid-IR. Our test demonstrates that we are seeing differences in 
the galaxy populations in both the global SED and the spectral in¬ 
dices. We, thus, conclude that our results about the ages of MQGs 
at z = 1.0 - 1.5 are robust and not an artefact linked to the SED- 
fitting degeneracies. 

The longer timescales for the senior galaxies, which are also 
relatively massive, is not directly consistent with a typical concep¬ 
tion in the downsizing scenario which states that the most mas¬ 
sive galaxies formed their stars early but also very rapidly. A 
short timescale for the formation of massive galaxies has also been 
claimed to be necessary to explain the or-element enhancement seen 
in early-type massive galaxies (such as ellipticals) in the nearby 


Universe (see, e.g., Thomas et al. 2005 and Renzini 2006). Our es¬ 
timations of the timescales for the most massive galaxies are not 
extremely long (typically 400 Myr) and would be rightly consis¬ 
tent with the values needed to match the chemical abundances (less 
than 1 Gyr; see Thomas et al. 1999; Worthey et al. 1992; Worthey 
1998 and references therein). In addition, these massive galaxies 
may be the product of mergers (or clumps) involving progenitors 
with shorter timescales, but with different ages, maybe linked to 
small offsets (of the order of tens or a few hundred Myr) in the 
ignition of the star formation. The SFHs obtained for the merged 
systems (i.e., our galaxies) would then present longer timescales, 
but do not violate any constrain linked to o-element enhancement. 
In Sect. 4.3, we further discuss the SFH of our galaxies as a func¬ 
tion of the mass and age and the importance of selection effects. 

The variety in the ages of quenched galaxies at z = 1.0 - 1.5 
is in agreement with previous works such as Bedregal et al. (2013) 
and Belli et al. (2015). Both papers present a wide range of stellar 
population properties for similar samples as ours. Bedregal et al. 
(2013), analysing the G102-G141 WFC3 grisms of 41 massive 
(log(M/M 0 ) > 10.65 M 0 ) quiescent galaxies at z ~ 1.5, with ex¬ 
ponentially declining SFH found that they had short star formation 
timescales (r < 100 Myr) and a wide distribution in stellar ages (1- 
4 Gyr). Belli et al. (2015) analysed Keck LRIS spectra of quiescent 
galaxies at 1.0 < z < 1.6 with log (M/M 0 ) > 10.6 M 0 and con¬ 
cluded that there are two different quenching routes. The youngest 
galaxies (t 0 ~ 1 Gyr) show a fast quenching (r ~ 100 Myr), while 
the older galaxies (up to 4 Gyr old) show slowly-declining SFHs 
(r > 200 Myr). Note, however, that our r values and theirs are not 
directly comparable, since we used different parameterizations for 
the SFH (i.e., the r parameter for an exponential is not exactly the 
same as the r value for a delayed exponential). 

4.2 Dissecting the UVJ diagram: distribution of stellar 
population properties 

In Figure 11, we plot the UVJ colour diagram, which was the start¬ 
ing point of our sample selection, to study how the UVJ colours 
do actually correlate with the derived galaxy properties. The main 
properties of the UVJ- and sSFR-selected sub-samples are shown 
in Table 5. We recall that we refer to galaxies in the quiescent region 
of the UVJ diagram and without IR detection as (/KZ-selected, 
while the sSFR-selected have sSFR < 0.2 Gyr -1 but are located 
outside the passive UVJ region (or inside but detected in the IR, 
see Sect. 2.2). 

Concerning ages, there is not a significant difference between 
those derived for the (/VT-selected and sSFR-selected samples. 
The average ages are (f M )=l.l*^ and 1.2g| Gyr for each sub¬ 
sample, respectively. However, the senior galaxies seem to have 
redder V — J colours, with average values (V-J)=1.46j'32, than the 
mature population, <V-J>= 1.31 }‘i|- 

The timescales are more clearly segregated in the UVJ plot. 
The r of the (/W-selected galaxies are shorter than those of 
the sSFR-selected galaxies: typical ( t ) values are 60^° Myr and 
250j“ Myr, respectively. Comparing our results with those in 
Belli et al. (2015), who also presented stellar population properties 
within the UVJ diagrams, we find a similar trend: they found that 
galaxies with r < 100 Myr were distributed in a narrow region par¬ 
allel to the diagonal line defined in the UVJ passive region. 

Belli et al. (2015) also found that galaxies with higher dust 
attenuations (A v ~ 1 mag) were outside the passive region. We 
find typical (A v )=0.5qj and (A v )=0.7q 4 mag for the (/VT-selected 
and sSFR-selected galaxies, respectively, in good agreement with 
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D4000 

Figure 10. Left panel: Stacks for the senior (thick red line) and mature (thick blue line) populations normalized in the Mgyy ( upper panel) and D4000 regions 
(lower panel). The average error of the stacks is ~ 10%. We used the dispersion as uncertainties (actual errors calculated by propagating the observed flux 
uncertainties were negligible). We also plot the normalized models with the average properties of each sub-population (thin red and blue lines). Right panel: 
Mguy versus D4000 plane, with the indices measured in the stacks for the mature, intermediate, and senior galaxies (blue, green and red diamonds with error 
bars). For comparison, we also show the evolutionary tracks for 2 models with different star formation timescales (a SSP and a delayed exponential model 
with t= 1000 Myr) and two different dust attenuations (A v =0.0 mag and A v =1.0 mag). The empty circles represent the expected index values for each track at 
different ages (given in Gyr in the plot). The three coloured lines show the tracks predicted by the models with the average properties of each galaxy population 
(given in the legend). The filled coloured circles mark the indices values at different ages (0.5, 1.0, 2.0, 3.0, 5.0 Gyr, from smaller to larger symbols). The 
indices measured in the stacks are compatible with the average properties derived from the SED-fitting for each sub-population. 


Parameter 

Sample 

Number 

Qi 

Median 

Q3 

z 

UVJ 

N=65 

1.05 

1.17 

1.25 


sSFR 

N=39 

1.03 

1.15 

1.24 

to [Gyi'l 

UVJ 


0.9 

1.1 

2.0 


sSFR 


1.0 

1.7 

2.5 

r [Myr] 

UVJ 


20 

60 

200 


sSFR 


160 

250 

400 

A v [mag] 

UVJ 


0.1 

0.5 

0.8 


sSFR 


0.4 

0.7 

1.1 

Z/Z G 

UVJ 


1.0 

1.0 

2.5 


sSFR 


0.4 

1.0 

2.5 

log(M/M e ) 

UVJ 


10.3 

10.5 

10.8 


sSFR 


10.3 

10.7 

10.8 

hi [Gyr] 

UVJ 


0.9 

1.1 

1.7 


sSFR 


0.8 

1.2 

1.7 

sSFR sed [Gyr -1 ] 

UVJ 


< 10 -4 

< 10 -4 

10 -3 


sSFR 


10 -2 

10 -2 

10 -1 

R [Gyr] 

UVJ 


0.7 

0.9 

1.2 


sSFR 


0.3 

0.6 

1.0 


Table 5. Galaxy properties of our sample according to their selection cri¬ 
teria {UVJ/ sSFR, see Sect. 2.2). The parameters are the same as those in 
Table 3. 


Belli et al. (2015). The percentage of (71/./-selected galaxies with 
A v > 1 mag is only 7%, while this percentage increases up to 28% 
for the sSFR-selected galaxies. All of the galaxies falling in the 
post-starburst region defined by Whitaker et al. (2011) (see Fig. 2) 


have A v < 0.5, i.e., they are relatively dust-free. Flowever, we do 
not see such a clear gradient in A v with the distance to the UVJ di¬ 
vision line as in Belli et al. (2015). A possible explanation for this 
discrepancy is the difference in the metallicity values used in each 
work. While Belli et al. (2015) used a very narrow range in metal¬ 
licity (a normal prior centered on the solar value Z=0.02 with a 
width of 0.005), we use 3 discrete values (Z=0.02, 0.05 and 0.008). 
Due to the dust attenuation-metallicity degeneracy, the dust attenu¬ 
ation distribution is obviously affected by the assumed metallicity 
values. 

With respect to the metallicity, no clear trends are present 
for any of the two selection criteria. The fraction of f/VT-selected 
galaxies with sub-solar, solar and super-solar metallicity is 20, 40 
and 40%, respectively. For the sSFR-selected, the percentages are 
25, 46 and 30%. We note that, if the metallicity was a fixed param¬ 
eter (e.g., solar metallicity), the galaxy properties would be better 
segregated within the UVJ diagram, as the best-fitting model would 
be constrained by only 3 parameters instead of 4. 

In Figure 1 1, we also present how SED-based sSFRs correlate 
with the position in the UVJ diagram. The most quiescent galaxies 
are located in a similar region to the galaxies with low r values. In 
fact, ~ 60% of the (/W-selected galaxies have sSFR SE d < 10 -5 
Gyr -1 , while this only happens for ~ 8% of the sSFR-selected 
galaxies. 

We have also derived the time since quenching, defined as 
the time since the galaxy became quiescent using our defini¬ 
tion from Sect. 2.2, i.e., how much time has passed since the 
galaxy had sSFR sed < 0.2 Gyr -1 . The t/V/-selected galaxies have 
been dead, on average, for almost 1 Gyr, (t q )=0.9\f Gyr, while 
the sSFR-selected galaxies have been more recently quenched, 
(f,>=0.6j° Gyr. 
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0.8 1.0 1.2 1.4 1.6 1.8 0.8 1.0 1.2 1.4 1.6 1.8 

V - J V - J 


Figure 11. UVJ diagrams with galaxies colour-coded on the basis of different properties, from left to right, top to bottom: mass-weighted age (?m), star 
formation timescale (t), dust attenuation (A v ), metallicity (Z), SED-based sSFR (sSFRsedX and time since quenching (t^). The circles represent galaxies with 
no degenerate solutions, while the triangles are galaxies with more than one solution in the Monte Carlo simulations (see Sect. 3.2). £/V/-selected galaxies 
are located inside the quiescent region of each diagram, while sSFR-selected galaxies fall outside by definition (see sect. 2.2, Table 5). 
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4.3 Tracing back the SFHs of MQGs: clues about their past 

In this Section we analyse the SFHs for MQGs at 1.0 < z < 1.5 
in the observed redshift range and at earlier lookback times. One of 
the advantages of determining SFHs with stellar population models 
is that they provide us with a time-dependent evolution of parame¬ 
ters such as the SFR and the galaxy stellar mass. Therefore, assum¬ 
ing closed-box evolution (i.e., no merging events or re-activation of 
the star formation activity), we do not only characterize the prop¬ 
erties of galaxies at the epoch of observation, but we can also trace 
back their properties and see how they assembled their mass at 
earlier cosmic times. At this point, we want to stress that the re¬ 
sults hereafter presented are strongly affected by the assumption of 
the SFH as a single burst, delayed exponentially declining model. 
This SFH parameterization can account for the star formation tak¬ 
ing place after one gas-rich major merger, but neglects the stellar 
population of the mass previously assembled (see Figure 1 from 
Hopkins et al. 2008). ACDM models and observed merger rates 
(e.g. Blucketal. 2009; Newman 2013; Man et al. 2014) predict 
that massive galaxies undergo at least one major merger since z=3. 
However, assuming multiple bursts of star formation, would imply 
the analysis of two stellar population models. This would signifi¬ 
cantly increase the degeneracies in the SED-fitting procedure and 
would severely complicate the interpretation of the results. We dis¬ 
cuss the impact in our results of a more complicated SFH in Ap¬ 
pendix A, but we delay a comprehensive analysis of more compli¬ 
cated SFH parameterizations to future papers. 


4.3.1 The SFR-Mass relation 

One of the most fundamental relations between galaxy properties 
is the SFR-Mass relation, the so-called MS. In Figure 12 we show 
the evolutionary paths of the galaxies in our sample in the SFRsed- 
Mass plane as a function of time. For this plot, we assumed that 
all galaxies started their formation at the same time (i.e., all SFHs 
share the t=0 point). The shape of the track depends mainly on the 
final mass and the star formation timescale, r: galaxies with short 
timescales form the bulk of their mass very quickly, they arrive 
relatively early to the MS and then their star formation decreases 
while their stellar mass remains almost unchanged (falling verti¬ 
cally). Galaxies with very long r values continuously increase their 
mass without changing their SFR (almost horizontal tracks) and 
may stay in the MS up to ~ 2 Gyr. According to our best-fitting 
SFHs (note that the choice of the parameterization is also impor¬ 
tant), after 100-500 Myr of evolution our galaxies would form a 
MS with a very similar slope to that observed directly at z > 1 with 
samples of star-forming galaxies (see, e.g., Speagle et al. 2014). We 
remark that, for this exercise, we offset the SFHs of our sample to 
make them match at t=0. The effect of galaxies starting their forma¬ 
tion at different epochs would result in a widening of the MS shown 
in Figure 12 for different ages. According to our results, galaxies 
would come out of the MS (lcr below) after approximately 1.5 Gyr 
(96% of our MQGs are located below the MS at that age), and 
quenching would then proceed almost vertically in this plot. We 
remark that this statement considers the main sequence evolution 
from z ~2 to z ~1. Therefore, we find that MQGs would pass a 
fair fraction of their lifetimes in the MS, with the possibility to live 
above the MS (in the starburst locus) for short periods of time of 
the order of 100 Myr. We conclude that the SFHs determined for 
the MQGs at 1.0 < z < 1.5 are consistent with the slope and even 
the location of the MS at z > 1 and that the existence of the MS for 


SFR SED 

> MS+2cr 

> MS+lcr 

MS ±lcr 

< MS-lcr 

< MS—2cr 

z bin 

i 

U) 

o 

N=2 

N=7 

N=8 

N=15 

N=12 


7% 

23% 

27% 

50% 

40% 

1.3-2.0 

N=8 

N=10 

N=5 

N=52 

N=46 


12% 

15% 

7% 

78% 

69% 

1.0- 1.5 

- 

- 

- 

N=104 

N=101 


- 

- 

- 

100% 

97% 

o 

oo 

1 

to 

- 

- 

- 

N=104 

N=104 


- 

- 

- 

100% 

100% 


Table 6. Number and percentage of galaxies above (1 or 2 <x), in and below 
(1 or 2<t) the MS from Speagle et al. (2014) at each redshift bin, as show in 
Figure 13. Note that the sum of percentages at each redshift does not equal 
100 as the ± hr bins also include the galaxies at ± 2<r. 

our sample of MQGs is mainly an age effect, in the sense that the 
MS is formed by galaxies with similar ages (~ 0.5-1 Gyr). 

We also want to remark that the current SFR sed for many of 
our galaxies derived from the SED modeling (averaged over the 
last 100 Myr of their history) is much smaller than the SFR ob¬ 
tained from observables. The current SFR S ed are n °t fully consis¬ 
tent with the SFR W and SFR //; measurements explained in Sec¬ 
tion 2.1 and used in the sample selection. Indeed, 45% of the sam¬ 
ple have SFR sed < 1O~ 3 M 0 yr~' but the lowest SFR estimated from 
the typical tracers is 10 0 5 M e yr _1 . This could be due to an over¬ 
estimation of the dust attenuation from the IRX-/? relation, but the 
choice of a delayed exponentially decreasing SFH is most proba¬ 
bly responsible for this effect. Typically, this parameterization pro¬ 
duces very low SFRs for dead galaxies (whose emission is dom¬ 
inated by an evolved stellar population), but there might exist a 
second stellar population with some (very low and negligible in 
terms of mass and emission) on-going star formation which would 
not be possible to recover with the assumed SFH parameterization. 
However, the fits of the model to our data are very good for all 
the sample and the inclusion of a second population would multi¬ 
ply by 2 the uncertainties and degeneracies. We, therefore, assume 
that our galaxies are dominated by the older stellar population and 
that one (composite) stellar population model explains the main 
features of our galaxies. We caution the reader, however, that the 
current SFRsed should be taken as lower limits. 

Following Kennicutt (1998), we can derive the SFR cor¬ 
responding to LIRGs/ULIRGs: L ;r = 10 11 L© corresponds to 
12 M 0 yr _1 , and L/ R =10 12 Lo to 121 M 0 yr~' (transformed 
to Kroupa 2001 IMF). All of the galaxies in our sample had 
SFR peaks larger than the LIRG limit (except one which has 
log(M/M 0 )= 10.0). The typical fraction of their lifetime spent in 
the LIRG phase is ~ 32% (~ 500 Myr). A significant fraction of the 
sample (46%) had SFR peaks larger than the ULIRG limit, but the 
typical fraction of their lifetime in this phase is much smaller, ~ 8% 
(~ 100 Myr). Our results favor LIRGs at z > 1.5 as the most likely 
progenitors of MQGs at 1.0 < z < 1.5, and that most of their stars 
were formed in star-forming events with SFRs around 100 M Q yr~*. 
The ULIRGs seem to be the progenitors of the most massive galax¬ 
ies. The fraction of galaxies which have undergone a ULIRG phase 
increases up to 65% (75%) when considering galaxies with masses 
larger than 10 10,8 M 0 (10 110 M 0 ). We discuss the connection be¬ 
tween MQGs and ULIRGs on the basis of their number densities in 
Sect. 4.3.2. 

In Figure 13, we show where our galaxies would be placed 
in the SFR SE D-Mass plane at different redshifts according to their 
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Figure 12. Evolutionary tracks in the SFRsED-Mass plane for our sample of MQGs. Left panel: the small coloured dots mark the location of galaxies at 
different times after their formation, as indicated in the legend. Note that these are direct ages from the SED-fitting and that we have assumed the same origin 
for all galaxies. We mark as thin grey lines the past path for each galaxy. The thick grey lines show the location of the MS at z ~ 1.2 (dashed-dotted line) and 
Z ~ 2 (solid line) according to Speagle et al. (2014), and at z ~ 1 according to Elbaz et al. (2007) (dashed line). The green shadow is the I cr dispersion of the 
MS at z=2.0. Most of the galaxies are ~ 0.5 - 1.0 Gyr old when they are on the MS region at z > 1.0. In fact, 75% of our sample is lcr below the MS at 
z=l .2 when they are 1 Gyr old and all of our galaxies are at least lcr below the MS, considering the MS evolution from z ~2 to z ~1, after 2 Gyr. Right panel: 
Location of our galaxies in the SFRsED-Mass plane at the epoch of observation (large filled circles), colour-coded by their best-fitting ages (to). We note that 
both the ages and SFRs used in this plot are the output of the SED-fitting (and not SFRuv neither tM )■ To better visualize the location of the galaxies in the 
SFRsED-Mass plane, we plot galaxies with very low SFRsed (< 10~ 3 M o yr 1 ) around SFRsed ~ 10 3 \f e yr 1 . 


past and future evolutionary tracks. We study 4 redshift bins cor¬ 
responding to 2 and 1 Gyr before the observations, the actual ob¬ 
served redshift and 1 Gyr after the observations. In Table 6, we 
show the percentages of galaxies above, within, and below the MS 
at different redshifts. Figure 13 shows that at 1.0 < z < 1.5, the 
epoch of observation, all of our galaxies lie lcr below the MS (this 
happens also when using the SFRs based on classical tracers), as 
expected given that we started with a selection of passive galax¬ 
ies. If we move 1 Gyr backwards in time, ~ 22% of the galax¬ 
ies are located in or above the MS, but a significant number frac¬ 
tion (~ 78%) of the galaxies at z ~ 1.6 are located lcr below the 
MS, meaning that they were already quiescent or in the process of 
quenching. At z ~ 2.1, the bulk of the galaxies cannot be plotted 
because ~ 70% have best-fitting ages smaller than 2 Gyr, so at that 
epoch their masses were much smaller than our limits in the plot 
(or they were not even formed yet). At that redshift bin, half of 
the galaxies (50%) are in the MS or above but the other half (15 
galaxies) are lcr below the MS as early as z ~ 2. In the opposite 
time direction, if we move 1 Gyr after the observations, all of our 
galaxies have very low SFRsed and are completely dead and well 
below the MS from Speagle et al. (2014) and Elbaz et al. (2007), as 
expected from purely passive evolution. 

Comparing with the parent population of galaxies more mas¬ 
sive than 10 10 M 0 at each redshift bin, we find fractions of quies¬ 
cent galaxies of 25, 12 and 3% at z = 1.0 - 1.5, z = 1.3 - 2.0 and 
z = 1.7 - 3.0, respectively. Different studies have shown that the 
fraction of quiescent galaxies constantly increases with time. For 
example, Dominguez Sanchez et al. (2011) derived a fraction of 
quiescent galaxies (log (M/M 0 ) > 10.6, sSFR < 0.01 Gyr -1 ) of ~ 20 
and 10% at z = 1.4-1.6 and z = 1.6-2.0 respectively. The percent¬ 
age found in Ilbert et al. (2013) for a sample of quiescent galaxies 
selected by their NUV-r vs r-J colours with log(M/M 0 ) > 9.6 in¬ 
creases from 6% to 13% from z = 2.5 - 3.0 to z = 1.1 - 1.5. The 
evolution of the fraction of quiescent galaxies (UVJ selected) in 
Muzzin et al. (2013) is not so strong (28% at z = 1.0- 1.5 and 24% 
at z = 2.0 - 2.5), although the mass limits used in this calculation 


change for each redshift bin, log(M/M 0 ) > 9.48 at z = 1.0 - 1.5, 
log(M/M 0 ) > 10.54 at z = 2.0 - 2.5. The fractions of MQGs that 
we derive at z > 1.5 by studying the past evolution of MQGs at 
z= 1.0-1.5 predicted from their SFHs are consistent with observa¬ 
tional results at higher redshifts. 


4.3.2 Number densities 

In the next paragraphs, we compare the number density of quies¬ 
cent galaxies in our work with a compilation of estimations from 
the literature. The number density of quiescent galaxies in our work 
for z > 1.5 are derived considering our results for the SFHs of 
1.0 < z < 1.5 MQGs and the position of the galaxies in a SFRsed- 
Mass plot at different epochs (Figure 13). In this case, we assume 
that galaxies evolve passively once their star formation is quenched 
and they do not experience any merger event. We warn the reader 
that a fair comparison of number densities requires taking into ac¬ 
count the differences in quiescent fraction linked to the stellar mass 
cut or the definition of quiescence, which vary from paper to paper 
and are difficult to consider in their full extent. Here we consider a 
galaxy as quiescent when is located lcr below the MS at each red¬ 
shift. With this method, the number densities of MQGs from our 
work are p=(7.0 ± 0.7)xl(T 4 Mpc -3 at z = 1.0 - 1.5, p=(2.3 ± 
0.3)xl0~ 4 Mpc- 3 z = 1.3 - 2.0, and p=(0.31 ± 0.08)xl0- 4 Mpc- 3 
z = 1.7-3.0. This is in good agreement with the number densities of 
quiescent galaxies reported in Muzzin et al. (2013) atz = 1.0 - 1.5 
(p=7.6xl0- 4 Mpc' 3 ), but smaller than their number densities at 
higher redshifts p=(3.3, 1.2, 0.65)xl0- 4 Mpc’ 3 at z = 1.5 - 2.0, 
z = 2.0 - 2.5 and z = 2.5 - 3.0, respectively. There are several fac¬ 
tors that could be affecting this comparison. In the first place, the 
redshift bins are not exactly the same. Besides, Muzzin et al. (2013) 
selection is only based on UVJ colours, while we are considering 
as quiescent galaxies lcr below the MS at each redshift, as derived 
from the backwards evolution of the MQGs at z = 1.0 - 1.5. We re¬ 
call that we have eliminated from our MQGs at z = 1.0 - 1.5 galax¬ 
ies detected in the IR in the quiescent UVJ region, when they have 
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Figure 13. Evolution of the MQGs in the SFRsed Muss plane at different redshifts. The observed properties (those measured directly from the data) are 
plotted as black filled circles in the redshift bin z = 1.0 - 1.5. The expected location of the galaxies 1 Gyr (z = 1.3 — 2.0) and 2 Gyr (z = 1.7 — 3.0) before the 
observations are plotted as black filled circles in the upper right and upper left panels, respectively. The lower right panel represents the expected location of 
the galaxies 1 Gyr after the observations (z = 0.8 - 1.2). We use the age, timescale, mass and SFR given by our SED fits to move forward and backwards in 
time in the SFRsED-Mass plane, assuming passive evolution. For comparison, we plot as a green line the MS from Speagle et al. (2014). The dark and light 
green areas mark the 1 and 2<r scatter, 0.2 and 0.4 dex. The blue line in the z = 0.8 — 1.2 panel is the MS from Elbaz et al. (2007). The orange/red lines show 
the SFR limit for FIRGs/UEIRGs converted into SFRs using Kennicutt (1998) relation. The light blue dotted line is the value of constant sSFR used in the 
sample selection (see Sect. 2.2; sSFR=0.2 Gyr -1 )- Galaxies with SFRsed < 1O -3 M 0 yr 1 are plotted as black empty diamonds at SFRsed ~ 10 1 M Q yr 1 . 


sSFR > 0.2 Gyr. We also note that galaxies in the quiescent region 
of the UVJ diagram at z = 1.7-3.0 can have much larger sSFR (up 
to 1.0 Gyr -1 ) than our quiescent limit. We recall that the values de¬ 
rived for the past evolution of MQGs must be taken with care, as we 
are assuming pure passive evolution and no merger events (which 
could help to decrease the number density of quiescent sources). 

We also derive the number densities of our observed galax¬ 
ies (z=1.0 - 1.5) for the different populations defined in Sect. 4.1 
(mature, intermediate, senior). This is shown in Figure 14, together 
with the number density of Muzzin et al. (2013) at z=1.0 - 1.5 
and Whitaker et al. (2013) at z=1.4 — 2.2. The existence of a pop¬ 
ulation of old (t > 1.3 Gyr) galaxies at z ~ 2 has already been 
found by Whitaker et al. (2013) after analysing the stellar popu¬ 
lations of 171 massive (/W-selected galaxies using stacked 3D- 
HST spectra. Given that in this work they use SSP models and 
we use more general (and realistic) delayed exponential models, 


the comparison with our sample is not straight forward, but we 
compare their results with our mass-weighted ages. Taking into 
account that there is a difference of ~ 1.6 Gyr between the me¬ 
dian redshift of Whitaker et al. (2013) sample, (z) = 1.64, and the 
median redshift of the senior population, ( z ) = 1.06, the galax¬ 
ies observed by Whitaker et al. (2013) would be ~ 3 Gyr old at 
the redshifts probed in our work, should they stay passive and 
not re-ignite. In our sample, we find a number density of galaxies 
with 1m > 2.0 of (1.1 ± 0.2) xlO -4 Mpc -3 for log(M/M 0 ) > 10.0 
and (0.9± 0.2) xlO -4 Mpc -3 for log(M/Me) > 10.5. These re¬ 
sults are consistent within errors with the number density of UVJ 
galaxies from Whitaker et al. (2013): (0.8 ± 0.1) xlO -4 Mpc -3 at 
z = 1.4 - 2.2 and with logfM/Mo) > 10.50. The similar number 
densities from Whitaker et al. (2013) for MQGs with ages around 
1.3 Gyr at z ~ 2 and those obtained in this work with ages > 2 Gyr 
at z ~ 1.1 is a good consistency check on the existence of old 
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★ This work, all ages 

★ This work, different ages 

■ Passive Muzzin+13 @ z= 1.0-1.5 
A Passive Whitaker+13 @ z= 1.4-2.2 


Figure 14. Number density of MQGs galaxies at z = 1.0-1.5 as a function 
of their mass-weighted age tM (red stars). The large black star is the total 
number density of MQGs from our work. We also show the number density 
of quiescent galaxies at the same redshift from Muzzin et al. (2013) (grey 
area). The small empty orange triangle depicts the number density and mean 
age of MQGs at z= 1.4 - 2.2, as measured by Whitaker et al. (2013). As¬ 
suming passive evolution for these galaxies, they would move to the orange 
filled triangle. The number density uncertainties (shown as 2<x) have been 
calculated assuming Poisson statistics, while the error bars in age represent 
the standard deviation in each age bin. 


galaxies at z > 1.5. However, this value is smaller than the number 
density of Muzzin et al. (2013) at z = 1.5 - 2.0 for galaxies with 
log(M/M 0 ) > 10.50 (p ~ 1.9xl0~ 4 Mpc~ 3 ). This difference could 
be due to the more restrictive selection from Whitaker et al. (2013), 
where the authors require G141 WFC3/HST grism detection for 
their sample. Given the discrepancies on previous results and the 
difficulty in making a fair comparison between number densities, 
it is difficult to reach firm conclusions regarding the assembly of 
the red sequence. However, the number densities (and fractions of 
quiescent galaxies) that we find when considering passive galaxies 
at z=1.0 - 1.5 and moving back in time seem to be consistent or 
smaller than those reported in studies based on galaxies at those 
actual redshifts above z=1.5. This could suggest that mergers play 
an important role or that the SFH of galaxies could be more com¬ 
plicated than the delayed exponentially declining assumed in this 
work (i.e., with the re-ignition of the star-formation for some galax¬ 
ies). 

To check the possibility that ULIRGs are the progenitors of 
MQGs, as mentioned in Section 4.3.1, we compare them in terms 
of number densities. Magnelli et al. (2013) found a number density 
p ~ 1.0xl0~ 4 Mpc~ 3 for ULIRGS at z = 2. Taking into account 
that 46% of our galaxies may have undergone a ULIRG phase, this 
would mean a number density of p ~ 3.2xl0~ 4 Mpc~ 3 . However, 
the duty cycle of the ULIRG phase is very short (— 100 Myr, see 
Sect. 4.3.1), which would significantly reduce the observed number 
density of ULIRGs at z=2. We conclude that the the possibility that 
ULIRGs z = 2 are the progenitors of the MQGs z = 1.2 is consistent 
in terms of number densities and timescales in rough terms. 


4.3.3 Age evolution 

In Figure 15 we show the evolution of the ages of MQGs over 
the last 10 Gyr. We compare our results with previous works 
based on stacked spectra (Mendel et al. 2015; Whitaker et al. 2013; 


Onodera et al. 2015; Schiavon et al. 2006; Choi et al. 2014) and in¬ 
dividual spectral measurements (Toft et al. 2012; Krogager et al. 
2013; van de Sande et al. 2013; Bedregal et al. 2013; Belli et al. 
2015; Barro et al. 2015). We compare the data with the predic¬ 
tions from pure passive evolution of a delayed r-models with dif¬ 
ferent formation redshifts. The median ages derived for our sample 
(Jm ~ 1.2 Gyr) are consistent with those front Belli et al. (2015) at 
the same redshift. The ages derived for the senior population could 
be explained in terms of passive evolution of the galaxies stud¬ 
ied at higher redshifts in Toft et al. (2012); Krogager et al. (2013); 
van de Sande et al. (2013); Whitaker et al. (2013); Onodera et al. 
(2015); Mendel et al. (2015) and Barro et al. (2015), suggesting 
that at least part of the quiescent population at z > 1.5 does not 
restart the star formation once they are quenched. However, the 
number density of old (age > 2 Gyr) MQGs is small and the av¬ 
erage properties of MQGs at z ~ 1.2 are dominated by galaxies 
with age < 2 Gyr. 

The observed properties of z < 1 galaxies are not fully consis¬ 
tent with a pure passive evolution of our sample of MQGs at z ~ 1.2. 
If the mature population evolved passively, they would have typi¬ 
cal ages consistent with those derived for the most massive galaxies 
(log (M/M 0 ) ~ 11.3) of the Choi et al. (2014) sample at z < 1. How¬ 
ever, the median mass of the mature population is ~ 10.4, meaning 
that they should grow by almost 0.9 dex in mass in 3-5 Gyr to 
be consistent with Choi et al. (2014) results. Recent works suggest 
that galaxies can increase their mass by a factor of ~ 2 due to resid¬ 
ual SFR (e.g. Perez-Gonzalez et al. 2008; Fumagalli et al. 2014) or 
via major mergers (Shankar et al. 2015), which is not enough to 
account for this mass discrepancy. The masses of the senior popu¬ 
lation (log((M)/M o )=10.7) would be consistent with the less mas¬ 
sive or intermediate mass sample of Choi et al. (2014), suggesting 
that they could be the progenitors of z < 1 MQGs. In this case, 
some level of rejuvenation should take place since z ~ 1.2, in the 
form of residual star formation or via mergers of senior and ma¬ 
ture (or even star forming) galaxies, in order to reconcile the ages 
of our sample with those from Choi et al. (2014). For example, 
~ 50% of the total mass should be formed in a burst of ~ 0.5 Gyr 
to make the ages of our senior galaxies (t M =2.6 at z=1.06) con¬ 
sistent with the ages observed by Choi et al. (2014) (t M =2.5 Gyr 
at z=0.6) for the intermediate mass population. The resulting to¬ 
tal mass would be log (M/M o )=11.0, also consistent with the val¬ 
ues observed by Choi et al. (2014). SFH with longer star formation 
episodes (t > 400 Myr) could also help alleviating the discrepan¬ 
cies. The properties of the less massive (log (M/Mq) ~ 10.7) quies¬ 
cent galaxies at z < 1 from Choi et al. (2014) cannot be explained 
in terms of purely passive evolution of quiescent galaxies at higher- 
z, suggesting that the less massive population is still forming stars 
at z > 1. However, the lack of number density considerations com¬ 
plicates the comparison with Choi et al. (2014), and may explain 
some of the discrepancies. 

The relatively uniform ages (around 1 -2 Gyr) measured at z > 1 
suggest that the quiescent galaxy population is being kept young by 
the constant addition of recently quenched galaxies, i.e., “new ar¬ 
rivals”, as we mentioned in Sect. 4.1 and was already stated in pre¬ 
vious studies (e.g. Mendel et al. 2015). We conclude that the forma¬ 
tion of the red sequence of quiescent galaxies is actually occurring 
at z = 1.0 - 2.0 (no results are available beyond that redshift). At 
these redshifts, the number density of the oldest population is small 
and the population of dead galaxies is dominated by new arrivals 
with ages around 1-2 Gyr (or mature galaxies, as defined in this 
work). Only at redshifts below z ~ L the MQG population is totally 
assembled and evolves passively with no significant new additions. 
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Figure 15. Evolution of the median stellar age of MQGs over the last ~ 10 
Gyr. The results from this work (red stars) are plotted for three l\j bins 
at their median redshift, the size of the symbols being proportional to the 
number density. The largest symbol is at tu ~ 1.2 Gyr, the median age of 
our sample. Our results are compared with other ages from the literature at 
different redshifts measured from stacks using SSP models (Mendel et al. 
2015; Whitaker et al. 2013; Onoderaetal. 2015; Schiavon et al. 2006) or 
from individual spectra using r-models (Toftetal. 2012; Krogager et al. 
2013; van de Sande et al. 2013; Bedregal et al. 2013; Belli et al. 2015). We 
also plot the mean luminosity-weighted age for 3 quiescent galaxies and the 
luminosity-weighted age for a galaxy in the process of quenching derived 
in Barro et al. (2015). The results from Choi et al. (2014) at lower redshifts 
are divided in 3 mass bins, log(M/M G ) ~ 10.7, 11.0, 11.3 from lighter to 
darker blue. The grey lines show the age evolution of models with SFHs 
of the form SFR(t) oc t/T 2 xe~'f T , with t = 100 Myr and different formation 
epochs (zf = 3.0,2.0,1.4) and r=400 Myr and z./ = 2.0. The yellow area 
indicates the redshift range studied in this work. 


This is in agreement with previous studies supporting a significant 
evolution of the number density of MQGs at lower redshifts: by a 
factor of ~ 2 from z ~ 1.5 down to z ~0.7, epoch at which the bulk of 
this population seems to be definitively assembled (see Daddi et al. 
2005; Eliche-Moral et al. 2010; Davidzon et al. 2013; Prieto et al. 
2013; Mendel et al. 2015; Prieto & Eliche-Moral 2015). 


4.3.4 Mass and age dependence of the SFH 

To study the dependence of the SFH with mass and age, we divide 
our sample in 3 mass bins in log(M/M G ) units: 10.0 - 10.5, 10.5 - 
10.8, and 10.8 - 11.5; and also in 3 bins in t M : < 1 Gyr, 1-2 Gyr, 
>2 Gyr (corresponding to the mature, intermediate and senior pop¬ 
ulations introduced in Sect. 4). We then construct the typical SFH 
for each sub-sample using the median r and SFR„ rav values. The 
SFHs are normalized so that, by integrating them, we recover the 
median mass of each subsample. We also derive the mean age of 
the Universe when the galaxies were formed (t u-form )• We plot the 
results of the SFRsed evolution for the three galaxy sub-samples 
divided by age and mass in Figure 16, and we give the obtained 
values in Table 7. 

If we study the SFH of our sample divided in age bins (left 
panel of Figure 16), we find that the senior galaxies must have 
been formed, on average, when the Universe was only 2 Gyr old, 
while the intermediate and mature galaxies are formed when the 
Universe was 3.4 and 4.1 Gyr old, respectively. The age of the Uni¬ 
verse when the galaxies had their maximum SFR is ~ 2.3, 3.6. 


4.2 Gyr with median SFR„ rav =60, 110, 230 M G yr~* and median 
r=400, 200, 60 Myr for the senior, intermediate and mature galax¬ 
ies, respectively. This is a direct implication of the galaxy proper¬ 
ties derived for each population and summarized in Table 3. Ac¬ 
cording to this scenario, the senior MQGs would have been formed 
in a relatively young Universe in longer and not very intense star 
formation episodes, while, on the other hand, the recently quenched 
MQGs at z = 1.0-1.5 must have been formed at later times in shorter 
and more intense bursts of star formation. We must notice that there 
are important selection effects in these results. The mature popula¬ 
tion must have been formed relatively quickly to be able to quench 
their star formation in less than 1 Gyr (their maximal age, by def¬ 
inition), meaning that they must present short r values to satisfy 
our quiescent selection criteria. Instead, the senior population does 
not present any selection bias, and could, in principle, have short or 
long timescales. 

With respect to the SFH divided in mass bins, the age of the 
Universe when they were formed was 3.6, 3.4 and 3.1 Gyr for the 
less massive, intermediate and most massive galaxies, respectively. 
The less massive galaxies are formed ~ 0.5 Gyr later than the most 
massive ones. This difference in time is larger than our age uncer¬ 
tainties (~ 0.2 Gyr), supporting that the most massive galaxies were 
formed first in time, in agreement with the downsizing scenario 
(Cowie et al. 1996). Recent cosmological simulations also predict 
that the most massive galaxies form most of the stars before and 
quench earlier with time (Zolotov et al. 2015). We derive that the 
peak of the SFR occurred at t u ~ 3.7, 3.6, 3.4 Gyr with median 
SFR„ rav =90, 120, 210 Mg yr~* and median r=80, 200, 250 Myr for 
the low, intermediate and high mass sub-samples, respectively. 

We compare these results with the scenario proposed in 
Thomas et al. (2005, 2010) for the formation of massive galaxies 
based on data in the nearby Universe (z 4, 0.06), which estab¬ 
lished that the most massive galaxies (log(M/M 0 ) ~ 12) formed 
2 Gyr after the Big Bang in relatively short formation timescales 
(r ~ 200 Myr), while lower mass galaxies assembled later and 
had longer star formation episodes (r > 1000 Myr). The galax¬ 
ies in our sample are formed ~ 2 Gyr earlier than the results from 
Thomas et al. (2010) suggest. For example, the age of the Uni¬ 
verse at the formation epoch of galaxies with log(M/M 0 ) ~ 11 pro¬ 
posed by Thomas et al. (2010) is ~5 Gyr (z ~ 1.2), while we find 
~3 Gyr (z ~ 2). However, a direct comparison with Thomas et al. 
(2010) is not straight forward. We recall that we are limiting our 
study to z=1.0-1.5, which restricts the latest formation epoch to 
t u-form ~ 6.0 Gyr, while Thomas et al. (2010) sample considers lo¬ 
cal galaxies which could have been formed at z < 1. The mass- 
dependence of the star formation timescale found by Thomas et al. 
(2010) is not seen in our sample. In contrast, the more massive 
galaxies in our sample have typically longer r than the less massive 
galaxies. Our r uncertainties are quite large (~ 20%) and the mass 
range analysed in this work is limited to 1.5 orders of magnitude, 
which may be reducing the mass-dependence of the star forma¬ 
tion timescales. But, more importantly, the median SFHs that we 
plot in Figure 16 are strongly affected by our selection criteria as 
we are only considering quiescent galaxies at 1.0 < z < 1.5. This 
biases our sample against less massive galaxies with long star for¬ 
mation episodes (long r) which were formed at z ~ 1.2 but are still 
forming stars at the redshift studied in this work. These galaxies 
did not have time to become quiescent at the redshift studied and 
therefore are not included in our sample. This biases the average 
r towards lower values. The number density of quiescent galaxies 
with log(M/M o )=10.0 - 10.5 is 10 times larger at z < 1 than at 
z = 1.0 - 1.5, meaning that this is an important selection effect. 
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redshift redshift 

4.0 3.0 2.0 1.5 3.0 2.0 1.5 



Figure 16. Average SFHs (SFR versus time) for our sample of MQGs divided in 3 age bins (left panel) and 3 mass bins (right panel). The yellow shaded area 
indicates the redshift studied in this work. 
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Table 7. Median SFH of MQGs at 1.0 < z < 1.5 divided in 3 !m bins and 3 
mass bins. We show the age of the Universe when the galaxies were formed 
(t u-form), the star formation timescale (r) and the maximal SFR (SFR„ MV ) 
for each subsample. 


5 SUMMARY AND CONCLUSIONS 

In this work, we have analysed the SFH of 104 MQGs 
(log(M/M 0 )> 10.0)atz = 1.0-1.5. They were selected on the basis 
of their location in the quiescent region of the UVJ diagram (and no 
IR detection) or imposing a cut in the specific star formation rate, 
sSFR < 0.2 Gyr 1 . We constructed the best possible SEDs by com¬ 
bining SHARDS spectro-photometric data, HST/WFC3 G102 and 
G141 grisms, and multi wavelength ancillary data form the Rain¬ 
bow database. The SEDs were compared to delayed exponentially 
declining stellar population models and a Monte Carlo algorithm 
was used to characterize in detail the uncertainties and the inherent 
degeneracies in this kind of study. 

Our main conclusions are: 

• The combined use of SHARDS and WFC3/HST grism data rep¬ 
resents a significant improvement in the study of MQGs at z=1.0- 
1.5. These data allowed to characterize and even break degenera¬ 
cies in some cases. For ~ 50% of our sample, we obtained a single 
solution in our analysis of the stellar population properties. 

• Two spectral features, the Mg vv and D4000 indices, were found 
to be very useful to disentangle between the SED-fitting degenera¬ 
cies in the other half of the sample. The spectral resolution of the 
SHARDS and grism data allowed us to break these degeneracies 


in ~ 30% of the cases where more than one model is consistent 
with the data. 

• The MQGs at z = 1.0 - 1.5 present a wide range of proper¬ 
ties, which are correlated between them. We divide our sample 
in 3 sub-populations based on their mass-weighted ages. Mature 
galaxies present ages t M < 1 Gyr, intermediate objects have t M = 1- 
2 Gyr, and senior galaxies are the oldest systems, t M > 2 Gyr. 
We find that mature galaxies have, on average, relatively short 
timescales, (r)=602g° Myr, smaller masses, log((M}/M o )=10.4j°3, 
and large dust attenuations, (A v )=0.8q 5 mag. The senior pop¬ 
ulation has longer star formation timescales, (t)=4003qq Myr, 
larger masses, log((M>/M o )=10.7] ( ! ) g and lower dust attenuations, 
(A v )=0.4g j mag. The intermediate population has transitional prop¬ 
erties: <7 m>= 1-4]'^ Gyr, (t)= 2003° 0 Myr, (A v >=0.5°j mag and 
l°g«M>/M 0 )= 10.5 Iq;* . 

• The global population of MQGs at z=1.0 — 1.5 is dominated by 
new arrivals; 85% of the sample is younger than 2 Gyr (in mass- 
weighted age). The progenitors of MQGs at z=1.0 - 1.5 started to 
form significant numbers of stars only after z ~ 2. 

• The existence of such different populations has been tested by 
the measurement of spectral indices ( Mguv and D4000) on the 
stacked data of each population. The indices values are consistent 
with those predicted by the tracks of the stellar population mod¬ 
els with the average properties of each population, thus confirming 
that the derived properties are real and not a consequence of the 
degeneracies present in the SED-fitting. 

• The UVJ colour-colour diagram segregates very well best-fitting 
properties such as r or sSFRsed- The t/VT-selected galaxies have 
short star formation timescales (r ~ 60 Myr) and low sSFR sed 
(< 10 -4 Gyr -1 ). By complementing the MQGs with galaxies lo¬ 
cated outside the UVJ passive region but with sSFRuv < 0.2Gyr~', 
we recover massive galaxies (log(M/M 0 ) ~ 10.7) with old popula¬ 
tions Um ~ 1.2 Gyr) and longer star formation timescales (r ~ 250 
Myr). 

• Analysing the derived SFHs, we studied the evolution of the SFR 
and mass as a function of time. We find that the tracks predicted by 
our SFHs are consistent with the slope and even the location of the 
MS of star-forming galaxies at z > 1.0. 

• According to the SFHs that we derive, all the MQGs of our sam¬ 
ple were LIRGs in the past, and about half of the sample went 
through an intense and short ULIRG phase. The typical time spent 
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in the LIRG/ULIRG phase was ~ 500/100 Myr (~ 32/8% of their 
lifetime). The fraction of galaxies which have undergone a ULIRG 
phase increases up to 75% when considering galaxies with masses 
larger than 10 110 M 0 , implying that the high-" ULIRGs may be the 
progenitors of the most massive quiescent galaxies. 

• The number densities and ages derived for the senior MQGs are 
consistent with the passive evolution of quiescent galaxies at higher 
redshifts (z ~ 2), meaning that at least part of the quiescent pop¬ 
ulation at z ~ 2 does not restart the star formation activity once 
quenched. However, at these redshifts, the number density of the 
oldest population is small and the population of dead galaxies is 
dominated by new arrivals with ages around 1-2 Gyr, i.e., the for¬ 
mation of the red sequence of quiescent galaxies is actually occur¬ 
ring at z= 1.0 - 2.0. Only at redshifts below z~ 1 the MQG popula¬ 
tion is totally assembled and evolves passively with no significant 
new additions. 

• The median SFHs of our galaxies divided in 3 age bins sug¬ 
gest that the senior MQGs at z= 1.0 - 1.5 have formed most of 
their stars in a relatively young Universe (ty ~ 2 Gyr) in long 
(r ~ 400 Myr) and not very intense star formation episodes 
(SFR,„ av ~ 60 M 0 yr _l ), while the mature MQGs must have been 
formed at later times (ty ~ 4 Gyr) in shorter (r ~ 60 Myr) and 
more intense (SFR„, nv ~ 230 M 0 yr~‘) bursts of star formation. 

• The median SFHs of our galaxies divided in 3 mass bins 
(log(M/M o )=10.0-10.5, 10.5-10.8, 10.8-11.5) suggest that the most 
massive galaxies at z= 1.0 - 1.5 were formed when the Universe 
was ~ 3 Gyr in intense (SFR m , u ~ 200 M G yr~') and relatively 
long (r ~ 250 Myr) star formation episodes, while the less mas¬ 
sive galaxies are formed 0.5 Gyr later in shorter (r ~ 80 Myr) and 
less intense (SFR,„ nv ~ 90 M Q yr~') star formation processes. 
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APPENDIX A: IMPACT OF THE ASSUMED STAR 
FORMATION HISTORY PARAMETERIZATION 

In the following Appendix we discuss the repercussion of our 
methodology in the presented results. One of the main limitations 



Figure Al. Mass-weighted stellar ages derived with r-models (used 
throughout the paper) versus the ages derived from SSP models. The black 
thick line is the one to one relation, while the dotted lines mark the separa¬ 
tion between the different populations according to their ages (see Section 
4.1). There is a good agreement between the two ages, with a median rel¬ 
ative error of 13% and no offset. Galaxies which are classified as mature, 
intermediate or senior with each of the two age values are highlighted in 
blue, green and red, respectively. 

of the work presented in this paper is the assumption of our SFH pa¬ 
rameterization in the form of a single-burst delayed exponentially 
declining function. Another possible issue is the effect of surviving 
degeneracies in our results (after considering the analysis based on 
indices and IR emission presented in Sect. 3.2.) 

Our results are not affected when the most significant solu¬ 
tions are chosen (the ones with the highest probability in the Monte 
Carlo simulations) instead of selecting the primary solutions with 
the procedure explained in Sect 3.3. In fact, only 12% of the pri¬ 
mary solutions are not the most significant ones. The largest differ¬ 
ences between the primary solutions and the most significant ones 
in terms of the average properties of the MQGs correspond to the 
senior population. The most significant solutions favor older me¬ 
dian ages (( t M ) ~ 2.9 Gyr versus 2.6 Gyr), longer star formation 
timescales ((r) ~ 500 Myr versus 400 Myr) and lower dust atten¬ 
uations (( A v ) ~ 0.3 mag versus 0.4 mag). For the mature and 
intermediate populations, the average properties remain unchanged 
when using the most significant solution. 

We have also tested the effect of using SSP models instead of 
delayed exponentially declining models as done in several papers 
in the literature (e.g. Schiavon et al. 2006; Whitaker et al. 2013; 
Onodera et al. 2015; Mendel et al. 2015). In Figure Al we show 
the comparison between the mass-weighted ages obtained with our 
fiducial delayed T-models and the ages obtained with SSP models. 
The degeneracies are not so significant for the SSP models (given 
that r is a fixed parameter) and therefore we compare the SSP ages 
coming from the most significant cluster of solutions. In general, 
there is a good agreement between the two ages, with a median 
relative difference between the two methods of 13% and no offset 
(median At=0.0). The fraction of mature, intermediate and senior 
galaxies are 39, 45 and 15%, respectively, very similar to the frac¬ 
tions considered throughout the paper (37, 48 and 15%). However, 
these are not exactly the same galaxies in each case. We highlight in 
Figure Al the galaxies that belong to each population determined 
with either of the two methods. The percentage of mature galaxies 
which are considered mature with each of the two methods is 80%, 
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Figure A2. Upper left: Stellar population models with the average properties of the mature population (dark blue) and the senior population (dark red) 
normalized to one solar mass. The difference in luminosity for the two models is evident (a factor of 0.6) and reflects the stronger emission of the younger 
stars. We also show two models which are the combination of the senior plus the mature population for different mass fractions (light blue represents 20-80% 
of the mass formed in the senior-mature burst, orange is 50-50% of the mass formed in each burst). The two composite models are normalized between 0.9 and 
3 pm rest-frame to account for the same mass as the mature model. The stacked photometric data for the mature population is plotted as white circles (black 
histogram for the grism data) and the grey shadow represents the lcr errors. Middle left: Flux ratio of the composite stellar population models with respect 
to the mature model at SHARDS resolution (R~50). Different colours represent different mass fractions, as shown in the legend. The SFHs are shown in the 
lower left with the same colour code. The dark and light grey shadowed areas represent 1 and 2cr average errors for the mature stack flux, respectively. At 
most, 20% (30%) of the mass could have been formed in a previous burst (3 Gyr before) and be hidden in our mature stack within lcr (2cr) errors. On the right, 
we do the same exercise, but now considering a composite stellar population combining the mature (dark blue) plus the intermediate (dark green) average 
properties. In this case, the intermediate model emission is comparable to the mature one, meaning that less than 5% of the mass could have been formed ~1.4 
Gyr ago in a relatively long burst (r ~ 200 Myr) and be hidden in our mature stack within lcr errors. 


72% for the intermediate, and 62% for the senior population. The 
dust attenuation distribution obtained with the SSP models resem¬ 
bles that obtained with r-models. Regarding metallicities, the SSP 
models favor sub-solar metallicity (~ 50% of the sample), while 
only 23% of the galaxies are fitted with sub-solar metallicities with 
the r-models. Although the ages, dust attenuations and fractions 
obtained with SSP or r-models are comparable, we want to remark 
that the goodness of the SED-fitting, based on visual inspection, is 
better when using r-models. 

Another important assumption to take into account is that we 
have only considered one star formation burst model in our anal¬ 
ysis. This may be a too simplified parameterization of the SFHs 
of real galaxies. A multiple-burst model is probably more realistic 
for massive galaxies. We have performed a test to check whether 
more complex SFHs could reproduce the SEDs of our sample. In 
this test, we have started with the stack for the mature population of 
massive quiescent galaxies at z = 1.0 — 1.5, best fitted by a 0.8 Gyr 


old stellar population with short timescale (50 Myr) and 0.8 mag 
of visual attenuation. We have then calculated how much stellar 
mass formed in a previous burst could be hidden in the form of an 
older stellar population. For these older stars, we have assumed a 
3 Gyr old burst and a timescale of 400 Myr. This means that we 
are checking whether one of our senior galaxies could have experi¬ 
enced a recent burst 1 Gyr ago overshining a significant amount of 
the older stars. The test points out that a 3 Gyr old burst accounting 
for 20% of the total stellar mass of our galaxies would contribute 
to the total SED with fluxes within the typical photometric error of 
the stack (8%). A larger fraction of mass formed in an older burst 
is incompatible with the stack for the mature population since we 
would be able to measure distinctive colours and absorption indices 
within our uncertainties. This is shown in the left panels of Figure 
A2. We have done the same exercise considering a burst taking 
place 1.4 Gyr earlier with r = 200 Myr, corresponding to the av¬ 
erage properties of the intermediate population. Fess than 5% of 
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the mass formed in such a burst could be hidden in the stack of 
our mature population within errors (right panels of Figure A2). 
Therefore, we conclude that the mature stack is dominated by the 
younger stellar population and that the contribution from an older 
population could account at most for 20% (30%) of the mass at lcr 
level (2cr). On the other hand, if we consider the stack for the se¬ 
nior population, the contribution of a young burst (0.8 Gyr) would 
become comparable to the older population for a fraction of mass 
of ~1%. This is due to the stronger emission of the younger stars, 
which shield the older stellar population for small mass fractions. 
This confirms that the contribution of a younger burst to the total 
mass of the senior population is negligible and that we are actually 
observing galaxies formed more than 2 Gyr ago. 

Analyzing our SEDs in terms of more complicated SFFls (in¬ 
cluding several or extended bursts for which delayed decreasing 
exponentials would be a poor fit) would imply significantly more 
complex degeneracies. This task is far beyond the scope of this pa¬ 
per, which we have based on the typical parameterization used in 
the literature (SSPs, exponentials) and a discussion about the ro¬ 
bustness of the derived ages and timescales. We delay the discus¬ 
sion about extended SFFls for future papers. 
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